MHD simulations of small ELMs at low triangularity in ASDEX Upgrade

A Cathey; M Hoelzl; G Harrer; M G Dunne; G T A Huijsmans; K Lackner; S J P Pamela; E Wolfrum; S Günter; the JOREK Team; the ASDEX Upgrade Team; the EUROfusion MST1 Team

doi:10.1088/1361-6587/ac5b4b

1. Introduction

The ITER experimental thermonuclear fusion reactor is predicted to operate in high-confinement mode (H-mode), which is characterised by the quasi-periodic excitation of type-I ELMs (edge-localised modes) [1]. In H-mode there are reduced levels of turbulent transport in the edge of the confined region, thus forming a narrow transport barrier that creates a 'pedestal' in the density and temperature profiles. Type-I ELMs are macroscopic magnetohydrodynamic (MHD) instabilities that are destabilised by a high pressure gradient and/or high toroidal current density (and current density gradient). In standard ELMy H-mode, the pedestal rises (thus increasing $\nabla p$ which simultaneously increases $j_\mathrm{tor}$ by the formation of the bootstrap current) until type-I ELMs are excited. The appearance of type-I ELMs causes the pedestal to crash in a timescale of ${10^{2}\sim10^{3}~\mathrm{\mu s}}$ , and it is followed by an inter-ELM phase that lasts $10^1 \sim10^3~\mathrm{ms}$ . The large transient heat loads associated with these ELMs must be avoided in future tokamaks in order to achieve an acceptable divertor lifetime [2, 3].

Naturally, ELM-free (e.g. QH-mode, I-mode, EDA H-mode) and ELM-mitigated (e.g. RMP mitigation, pellet pacing, grassy ELMs, small ELMs) operational conditions have been successfully achieved using several methods in different existing tokamaks. However, naturally ELM-free regimes and ELM-mitigated scenarios are only achievable in reduced parameter regimes that differ between different tokamaks and, therefore, it is uncertain which methods and regimes will be accessible in future tokamaks [4]. Extrapolating such regimes to ITER becomes increasingly uncertain, because ITER parameters cannot be simultaneously accessed by existing tokamaks (e.g. density or collisionality, but not both at the same time) [5].

Certain small ELM regimes are completely free of type-I ELMs, maintain other desirable features, and could signify a fitting alternative for ITER. One such regime is the quasi-continuous exhaust (QCE) regime, which maintains good confinement properties (previously also called type-II ELMs [6, 7] or simply small ELMs [8, 9]), and is routinely operated in the ASDEX Upgrade (AUG) and TCV. In order to access the QCE regime, it is necessary to operate at high separatrix density ( ${n_{e,\mathrm{sep}}\gtrsim 0.3\times n_\mathrm{GW}}$ , where ${n_\mathrm{GW}}$ is the Greenwald density), with high triangularity and close to double null [9]. The physical mechanism that constrains the pedestal beneath the type-I ELM stability boundary remains unclear, but it is thought that high-n ballooning modes (where n is the toroidal mode number) that are located near the separatrix and cause sufficient outward transport of heat and particles, which impinge onto the divertor [8, 9]. It has been found that the power fall-off length is larger during the QCE regime than expected from the empirical Eich scaling [10]. Other small ELM regimes that can be free of type-I ELMs include grassy ELMs (observed in JT-60U [11], JET [12], EAST, and it is foreseen as a potential operational regime for CFETR [13]), type-III ELMs (are not reactor-relevant because they cause confinement degradation [14]), and a low density small ELM regime in JET [15]. Another type-I ELM-free regime that is now routinely operated in AUG is the enhanced D-alpha (EDA) H-mode [16] (first found in Alcator C-mod [17]). The first experiments of EDA H-mode in AUG were performed at low triangularity, with pure electron heating and observed a narrow operational window in terms of the applied heating power [16]. Recently, it has been found that the operational window can be extended to higher heating powers by increasing the plasma triangularity; heating above said operational window results in a type-I ELMy H-mode [18]. This operational regime always features an edge quasi-coherent mode (QCM) (in a frequency range between 20 and $80~\mathrm{kHz}$ ), which is thought to be responsible for regulating the pedestal below the type-I ELM stability boundary [19].

The JOREK non-linear extended MHD code [20, 21] has been extensively used to simulate macroscopic edge instabilities in tokamaks plasmas. In particular, for AUG, it has been used to produce realistic simulations of type-I ELMs [22, 23], RMP-ELM mitigation and suppression [24], and pellet-triggered ELMs [25, 26]. The present article details JOREK simulations of small ELMs at low triangularity in AUG and discusses their relation to the small ELMs that underlie the QCE regime and to the QCM of EDA H-mode. The simulations presented here follow on from an approach for modelling the pedestal build-up described in [22]. In simulations at sufficiently high separatrix density ( ${n_{e,\mathrm{sep}}\approx0.4\times n_{GW}}$ ), small ELMs appear beneath the type-I ELM stability boundary and feature resistive peeling-ballooning modes near the separatrix, which cause quasi-continuous heat exhaust. The present article is structured as follows. A brief description of different types of small ELMs in AUG and some features of the EDA H-mode are presented in section 2. The JOREK model used for the present simulations, together with the simulation set-up details and the axisymmetric pedestal build-up, are presented in section 3. In section 4, the results of the non-axisymmetric simulations are presented, and a detailed analysis is provided. Two different paths to leave the small ELM regime and reach a type-I ELMy H-mode are presented in section 5. Finally, conclusions and an outlook for future work are discussed in section 6.

2. Small ELMs and EDA H-mode at ASDEX Upgrade

The term 'small ELMs' has been used in AUG as a broad category that includes small amplitude ELMs (relative ELM size ${\lesssim} \ 5\%$ of the plasma stored energy), but excludes type-III ELMs. Type-III ELMs appear as distinct peaks in the divertor shunt currents with a roughly constant repetition frequency ( ${f_\mathrm{type-III}\sim10^3~\mathrm{Hz}}$ ) as opposed to the quasi-continuous activity of the QCE regime and EDA H-mode.

As mentioned in the introduction, the important ingredients for maintaining the QCE regime are high separatrix density ( ${n_{e,\mathrm{sep}}/n_\mathrm{GW}\gtrsim0.3}$ ) and closeness to double null (together with high triangularity) [6, 8, 9]. Given that there is little variation of separatrix temperature in a given device ( ${T_{e,\mathrm{sep}}\approx100~\mathrm{eV}}$ for H-modes in AUG [27]), high $n_{e,\mathrm{sep}}$ translates to high separatrix collisionality ( ${\nu_{e,\mathrm{sep}}^*\propto n_e/T_e^{\,2}}$ ). In existing tokamaks, high ${\nu^*_{e,\mathrm{sep}}}$ implies also a high pedestal collisionality, because the temperature cannot increase arbitrarily, due to the excitation of type-I ELMs. On the other hand, ITER can reach higher pedestal top temperatures and is expected to operate with high ${\nu^*_{e,\mathrm{sep}}}$ and low ${\nu^*_{e,\mathrm{ped}}}$ . Such conditions cannot be achieved simultaneously in existing tokamaks; therefore, fundamental uncertainties exist regarding whether ITER could operate in the QCE regime [8].

Depending on the triangularity and the edge safety factor, increasing heating power can either lead to a sustained small ELM regime (with 'standard' QCE parameters: high triangularity, high q₉₅, and close to double null) or to a transition from small ELMs to type-I ELMs (with 'standard' QCE parameters, but lower q₉₅), as shown in figure 1. For instance, discharge #35572 ( ${-2.0~\mathrm{T}}$ , ${1.0~\mathrm{MA}}$ , ${\delta = 0.397}$ , low safety factor ${q_{95} = 3.73}$ , ${P_\mathrm{heat} = 10~\mathrm{MW}}$ ) heated with ICRH and NBI ( ${P_\mathrm{rad}\approx4.1~\mathrm{MW}}$ ) shows small ELMs, and which transitions to a type-I ELMy H-mode upon increasing the heating power to ${P_\mathrm{heat} = 14.6~\mathrm{MW}}$ ( ${P_\mathrm{rad}\approx5.2~\mathrm{MW}}$ ). Conversely, at a higher q₉₅, discharge #39565 ( ${-2.5~\mathrm{T}}$ , ${0.8~\mathrm{MA}}$ , ${\delta = 0.401}$ , ${q_{95} = 5.76}$ ) heated with ECRH and NBI, retains the pure small ELM behaviour through a stepped increase in the heating power up to ${P_\mathrm{heat} = 13~\mathrm{MW}}$ ( ${P_\mathrm{rad}\approx5.5~\mathrm{MW}}$ ). Discharge #35572 is chosen as a comparison to #39565, despite their differences in heating power and plasma stored energy, because so far in AUG a transition from small ELMs to type-I ELMs, achieved by increasing the heating power, has not been obtained in the standard QCE regime shape. A quantitative experimental investigation of the heating power needed to transition the QCE regime to a type-I ELMy H-mode goes beyond the scope of this work.

**Figure 1.** From top to bottom, the heating power (radiated power in dotted lines), plasma-stored energy, and ELM monitor of two AUG discharges with small ELMs that feature increasing heating power. A discharge with low edge safety factor (left; #35572, ${q_{95} = 3.73}$ ) changes from a small ELM-dominant regime to a type-I ELM-dominant regime upon increasing the heating power. The discharge with high ${q_{95}}$ (right; #39565, ${q_{95} = 5.76}$ ), on the other hand, remains in a small ELM-dominant regime throughout the heating power steps.
Download figure:
Standard image High-resolution image

**Figure 1.** From top to bottom, the heating power (radiated power in dotted lines), plasma-stored energy, and ELM monitor of two AUG discharges with small ELMs that feature increasing heating power. A discharge with low edge safety factor (left; #35572, ${q_{95} = 3.73}$ ) changes from a small ELM-dominant regime to a type-I ELM-dominant regime upon increasing the heating power. The discharge with high ${q_{95}}$ (right; #39565, ${q_{95} = 5.76}$ ), on the other hand, remains in a small ELM-dominant regime throughout the heating power steps.
Download figure:
Standard image High-resolution image

Even with a single null configuration and low triangularity, small ELMs can be achieved at sufficiently high separatrix density. Nevertheless, such small ELMs are associated with degrading confinement, and even to an H-mode density limit (HDL) ( $n_e\lesssim n_\mathrm{GW}$ ). In such cases, as the separatrix density is increased (by increasing the gas puff rate), filamentary transport also increases [28, 29] and can lead to a flattening of the pressure gradient. This, in turn, causes a reduction of the edge radial electric field⁹ , which is what causes the back transition to L-mode, i.e. the HDL [28]. Recently, a correlation has been observed between a ballooning stability at the separatrix and the onset of the HDL in AUG [31] and in JET-ILW [32]. In this context, the deterioration of the H-mode confinement and, ultimately, the breakdown of the H-mode, are explained by an excess of cross-field transport caused by small ELMs located near the separatrix.

The enhanced D_α H-mode is a no-ELM operational scenario with high density, which is potentially attractive for ITER. Type-I ELMs are not destabilised during EDA H-mode operation because the pedestal is not able to build-up sufficiently. The transport mechanism that allows the pedestal to remain below the type-I ELM stability boundary is thought to be an electromagnetic mode dubbed QCM [19]. At low triangularity, the EDA H-mode lives in a very narrow operational space in terms of the applied heating power; however, in AUG this operational space has been observed to expand when increasing the plasma triangularity [18]. Nonetheless, a stationary EDA H-mode inevitably transitions to a type-I ELMy H-mode upon a sufficient increase in the heating power. In AUG, the QCM moves in the electron diamagnetic direction and has a fluctuation frequency in the range $f\approx20{-}80~\mathrm{kHz}$ [16]. Similarly, type-II ELMs are accompanied by broadband fluctuations (in magnetic pick-up coils and ECE signals) with frequencies in the range of $f\approx30{-}50~\mathrm{kHz}$ [7].

3. JOREK simulation set-up and axisymmetric build-up

The simulations presented in this paper were produced with the reduced MHD version of the 3D non-linear extended MHD code JOREK [20, 21]. JOREK is a versatile code with an implicit time-stepping scheme that can consider realistic X-point geometries by making use of a 2D finite element representation of the R-Z plane, together with a Fourier representation of the toroidal direction. The reduced MHD model in JOREK simplifies the visco-resistive MHD equations with two considerations. First, the toroidal magnetic field is constrained to be time-independent ( ${\boldsymbol B_\varphi = \frac{F_0}{R}\boldsymbol{\hat{\varphi}}}$ , where F₀ is a constant, R is the major radius, and $\boldsymbol{\hat{\varphi}}$ is the toroidal coordinate). Second, the poloidal velocity is considered to be comprised of the $\mathrm{ExB}$ velocity, which is further considered to lie in the poloidal plane, i.e. ${\boldsymbol v_\mathrm{pol} = \boldsymbol v_\mathrm{ExB}}$ ; this assumption allows a potential formulation for the poloidal velocity through the electrostatic potential ${\boldsymbol E = -\boldsymbol\nabla \Phi}$ . It is then possible to include diamagnetic effects (as an extension to the reduced MHD model) by considering the poloidal velocity to be

$\begin{align*} \boldsymbol v_\mathrm{pol} & = \boldsymbol v_\mathrm{ExB} + \boldsymbol v_{*,i} \\ & = -\frac{R}{F_0}\boldsymbol\nabla\Phi\times\boldsymbol{\hat{\varphi}} - \frac{m_i R}{e F_0~\rho}\boldsymbol\nabla p_i\times\boldsymbol{\hat{\varphi}}, \end{align*}$

where ${m_i}$ and ${p_i}$ are the ion mass and pressure respectively, e is the fundamental electric charge and ρ is the mass density.

The resulting model is a closed system of five equations for the poloidal magnetic flux (ψ), the electrostatic potential (Φ), the parallel velocity ( $v_\parallel$ ), the mass density (ρ), and the single fluid pressure ( $\,p = p_e+p_i$ ). Including diamagnetic effects allows the JOREK simulations to recover realistic radial electric fields in the pedestal region, i.e. ${E_r \propto n_i^{-1}\nabla p_i}$ [30], which play a fundamental role in the stability of PB modes (particularly those with high toroidal mode numbers) [33]. The stability of PB modes is also jointly determined by the edge current density (and its gradient), which is comprised of an Ohmic contribution and a bootstrap current contribution. These are included in JOREK by assigning a current density source determined by the initial current density profile (which is in turn comprised of Ohmic and bootstrap current contributions, ${j_0 = j_{0,\Omega}+j_{0,\mathrm{bs}}}$ ). The time-evolving contribution from the bootstrap current density is considered in JOREK by making use of the Sauter analytical expression [34, 35]. The total current density source then corresponds to ${j_\mathrm{source}(t) = j_0 + [\,\,j_\mathrm{bs}(t) -j_{0,\mathrm{bs}} ] = j_{0,\Omega} + j_\mathrm{bs}(t)}$ , and it is included in the induction equation

$\begin{align*} \frac{1}{R}\partial_t\psi & = - \eta (\,j-j_\mathrm{source}) - \frac{1}{F_0} [\Phi,\psi] -\frac{1}{R}\partial_\varphi\Phi \\ &\quad + \frac{m_i}{e F_0~\rho}\left( [\,p_e,\psi] + \frac{F_0}{R} \partial_\varphi p_e \right), \end{align*}$

where ${[A,B] = (\partial_R A)(\partial_Z B) - (\partial_Z A)(\partial_R B)}$ is the Poisson bracket, and η is the resistivity. The extension to include the bootstrap current density as a source term in JOREK was introduced in [36], and it has been used in recent simulations of ELMs [22, 25, 26]. The vectorial momentum equation of the extended MHD model,

$\begin{align*} \rho\partial_t \boldsymbol v & = -\rho\boldsymbol v \cdot\nabla\boldsymbol v -\rho\boldsymbol v_i^*\cdot\nabla\boldsymbol v_\mathrm{ExB} - \nabla p \\ & \quad + \boldsymbol j\times\boldsymbol B + \nabla\cdot\underline{\boldsymbol \tau} + \boldsymbol S_{\boldsymbol v}, \end{align*}$

where $\boldsymbol v = \boldsymbol v_\parallel + \boldsymbol v_\mathrm{ExB}$ , is solved by separating it into two equations through different projections; one of them by projecting it in the direction parallel to $\boldsymbol B$ and the other by projecting it in the poloidal plane. $\boldsymbol S_{\boldsymbol v}$ is a source term (equation (15) in [20]) and $\underline{\boldsymbol \tau}$ is the viscous stress tensor (equation (9) in [20]). Finally, the equations of the mass density and the single fluid pressure are given by equations (32) and (33), respectively in [20].

All of the simulations presented in this paper consider diamagnetic effects and the bootstrap current density source unless specified otherwise. The remainder of this section describes the numerical set-up of the simulations and the axisymmetric build-up of the pedestal profiles.

3.1. Numerical parameters and simulation set-up

Initial conditions that are stable to ideal peeling-ballooning modes are considered in this work. These are obtained from a post-ELM equilibrium reconstruction of AUG discharge #33616 at roughly ${7~\mathrm{s}}$ with the CLISTE code [37]. The magnetic field at the magnetic axis is $2.5~\mathrm{T}$ and the plasma pressure is ${I_p = 0.8~\mathrm{MA}}$ . We consider a lower single null magnetic configuration with low triangularity ${\delta_\mathrm{av} = 0.29}$ and with the ion ${\boldsymbol B\times\boldsymbol\nabla B}$ drift direction pointing towards the active X-point. The reconstructed current density profile is constrained based on the steepness of the density, temperature, and pressure profiles, i.e. with a bootstrap current constraint [38].

The initial profiles at the outboard midplane of the electron density, single fluid temperature and pressure, and toroidal current density are shown in figure 2. The plasma temperature, T, is the sum of the ion and electron temperatures, which are assumed to be equal in the physical model used for the simulations presented here. The toroidal current density is comprised of an Ohmic contribution together with a Pfirsch–Schlüter and a bootstrap current contribution¹⁰ . The bootstrap current has only a small contribution in the initial profiles, due to the small steepness of the initial pressure profile. An ideal MHD stability analysis with the MISHKA code [39] indicates that these post-ELM pedestal profiles are stable (to ideal PB modes) and would first need to steepen in order to excite a type-I ELM.

**Figure 2.** Initial conditions for the electron density and plasma temperature ( ${T = T_e+T_i}$ ) (a), and plasma pressure ( ${\,p = p_e+p_i}$ ) and toroidal current density (b) in the outboard midplane.
Download figure:
Standard image High-resolution image

**Figure 2.** Initial conditions for the electron density and plasma temperature ( ${T = T_e+T_i}$ ) (a), and plasma pressure ( ${\,p = p_e+p_i}$ ) and toroidal current density (b) in the outboard midplane.
Download figure:
Standard image High-resolution image

We use a flux-surface aligned grid that considers the confined region, the scrape-off layer, and the private flux region. Details of the methods used to construct flux-surface aligned grids in JOREK can be found in section 2.2 of [20]. The grid is comprised of 138 points in the radial direction (120 points in the confined region and 18 in the scrape-off layer) and 354 in the poloidal direction. A convergence scan in the grid resolution has been performed; we observe that the linear growth rates of instabilities with ${n\leqslant20}$ do not change by changing the radial and poloidal resolution (higher mode numbers were not probed because the dominant modes are ${n\leqslant12}$ due to diamagnetic stabilisation). For the axisymmetric build-up, only one poloidal plane (i.e. the R-Z plane at fixed ϕ) is considered, and for the non-axisymmetric simulations 32 poloidal planes are used in order to simulate the toroidal mode numbers ${n = 0,2,4,\dots,12}$ through a Fourier decomposition (and 64 planes for ${n = 0,2,4,\dots,20}$ ).

The resistivity at the magnetic axis that is used for the simulations is ${\eta_\mathrm{axis} = 6.6\times10^{-8}~{\Omega \textrm{m}}}$ , and it follows the Spitzer temperature dependency, ${\eta\propto T^{\,-3/2}}$ . This value is larger than the Spitzer resistivity along the magnetic axis ( ${\eta_\mathrm{Spitzer,axis}\approx2.1\times10^{-9}~{\Omega \textrm{m}}}$ ). This temperature dependence is incomplete in the pedestal due to neoclassical effects and impurities. The true experimental pedestal resistivity is larger than the Spitzer value, due to neoclassical effects and $Z_\mathrm{eff}$ (which increases from the core to the edge) such that the resistivity used in the simulations agrees with the experimental values within the error bars in the pedestal region, which is where the instabilities of interest are located. The parallel electron heat diffusion may be estimated with the Spitzer–Härm expression, ${\chi_{\parallel,\mathrm{SH}^e} = 3.6\times10^{29} T_{e,[\mathrm{keV}]}^{\,5/2}/n_{e,[\mathrm{m^{-3}}]}}$ [40]. For a plasma temperature of ${1~\mathrm{keV}}$ and density of ${5\times10^{19}~\mathrm{m^{-3}}}$ , i.e. ${\psi_\mathrm{N}\approx0.8}$ in the post-ELM equilibrium shown in figure 2, ${\chi_{\parallel,\mathrm{SH}^e}\approx 1.27\times10^{9}~\mathrm{m^2}\,\mathrm{s}^{-1}}$ . For the simulations presented in this paper, the parallel heat diffusion is roughly 15 times lower than the Spitzer–Härm value, as motivated by the heat-flux limit which can account for a reduction of ${\chi_\mathrm{SH}}$ by a factor of 15–130 [41]. Fast parallel heat transport acts as a stabilising agent for PB modes [42].

3.2. Axisymmetric pedestal build-up

Together with the post-ELM equilibrium reconstruction of AUG discharge #33616, a pre-ELM reconstruction from the same discharge was considered. An ideal MHD stability analysis (using the MISHKA code) shows that the pre-ELM profiles are unstable to ideal PB modes, and that the post-ELM profiles are stable. Starting from the post-ELM profiles, we impose ad-hoc diffusion coefficients (that describe a well in the pedestal region) and sources that drive the pedestal towards the pre-ELM profiles. The diffusion coefficients and sources are constant in time as the simulation progresses. This results in a pedestal build-up at fixed pedestal width, which is a simplification of what is experimentally observed. Time-evolving diffusion coefficients and sources would require including several key physical effects that are beyond the scope of MHD, and which will be investigated in future work (more on this subject in section 3.3). Despite such limitations, this approach at modelling the pedestal build-up at fixed pedestal width was already successful at reproducing type-I ELM cycles for the first time [22].

As the pedestal n_e and T evolve due to the stationary diffusion and sources, the radial electric field and j_ϕ become driven by the increasing influence of diamagnetic effects and the bootstrap current density, respectively. Figures 3(a)–(d) show the time evolution of the pedestal in terms of the electron density, plasma temperature, radial electric field, and toroidal current density. The colours of the profiles change gradually from purple to blue with increasing time as shown in the colour bar on top of the figure. The profiles are plotted every ${0.2~\mathrm{ms}}$ during the first ${4~\mathrm{ms}}$ of an axisymmetric simulation. In the first millisecond, the profiles change shape quickly (as can be most clearly evidenced in the evolution of the density and radial electric field). Section 4 will show that non-axisymmetric instabilities, driven by the steepening profiles, prevent the pedestal from building up towards the profiles shown in figure 3, and that the underlying instabilities are of a resistive nature. From linear stability analysis, we know that the pedestal build-up shown in figure 3 does not cross the ideal peeling-ballooning boundary in the time shown; several simulations with the linear ideal MHD stability code MISHKA [39], which considers only the confined region but excludes the separatrix, were performed at different time points during the pedestal build-up ( ${t = 0.0,~0.5,~1.0,}~\mathrm{and}~5.0~\mathrm{ms}$ ) and it was found that all time points are stable to ideal PB modes.

**Figure 3.** Outboard midplane pedestal profiles of electron density (a), plasma temperature (b), radial electric field (c), and toroidal current density (d) during the imposed pedestal build-up. The first millisecond is defined by a strong steepening of the pedestal, and the pedestal top grows progressively at a fixed width.
Download figure:
Standard image High-resolution image

3.3. Limitations of the present approach

The pedestal build-up considered for the simulations presented in this work considers fixed diffusion coefficients and, as such, it implicitly assumes a constant level of turbulent and neoclassical transport. Experimentally, however, turbulent and neoclassical transport in the pedestal is known to evolve in sub-millisecond and millisecond timescales, respectively. These changes in transport levels would determine how exactly the pedestal top and width evolve. In this sense, to produce a fully realistic pedestal build-up, which would include the pedestal widening, it is necessary to run neoclassical and gyro-kinetic (or even kinetic) simulations to determine the dynamical turbulent and neoclassical transport throughout the simulation time. Unfortunately, state-of-the-art gyrokinetic simulations are presently unable to routinely model pedestal dynamics in an affordable fashion and, as such, attempting to produce a fully realistic pedestal build-up model would represent not only an extremely costly endeavour from the point of view of computation time, but also would require significant efforts in terms of code development, which lie well beyond the scope of the present work.

The ion and electron species can be approximated to have the same temperatures and densities in the pedestal region, but in reality ions and electrons behave in distinct ways, due to the large difference in their respective masses. Such effects are neglected in the present simulations, since we use the single fluid version of the JOREK code. However, a two temperature model has been developed in JOREK and it will be used in the future to understand the effect of such temperature separation in ELM physics, and advance the level of realism in our simulations. Another important physical effect that is not considered in the present approach is the penetration of neutral particles into the confined region and their interaction with the plasma (more generally speaking, a more complete SOL/divertor model is missing). The ensuing ionisation of the neutral particles would determine the amount of particle fuelling that should be considered at any given time during a simulation. These fuelling effects, in turn, directly influence the density (and ultimately temperature) profiles in the pedestal. Ongoing efforts are underway that permit JOREK simulations to consider such effects either by a kinetic treatment [43] or a fluid treatment [44, 45] of the neutrals.

4. Non-axisymmetric simulations

This section details the simulation results by first describing the linear growth phase of the instabilities, together with the early non-linear phase (phase during which n ≠ 0 modes interact with each other but not with the n = 0 axisymmetric background), which takes place roughly during the first millisecond of simulation time. Then, in section 4.2 a scan of resistivity during the linear growth phase is presented. Resistivity is observed to have a strong impact upon the linear growth rates, which allows us to identify the underlying instabilities as resistive peeling-ballooning modes located near the separatrix. Apart from the results of section 4.2, all simulations presented in this paper consider a resistivity of ${\eta_\mathrm{axis} = 6.6\times10^{-8}~{\Omega \textrm{m}}}$ in the magnetic axis, with a Spitzer temperature dependency. Later on, in section 4.3, it is shown that it is critical to include diamagnetic effects in order to properly model the instabilities under consideration. In the following sections (4.4–4.6) the non-linear phase of the simulations is described in detail.

4.1. Linear and early non-linear phases

As mentioned before, the pedestal does not cross the ideal PB stability boundary during the axisymmetric build-up shown in figure 3. However, non-ideal instabilities can become excited due to the finite resistivity used in JOREK. A simulation that includes the toroidal mode numbers ${n = 0,2,4,\dots,12}$ is performed by introducing perturbations with said mode numbers at noise-level amplitudes after ${t\approx0.1~\mathrm{ms}}$ of axisymmetric build-up. The odd mode numbers were neglected for these simulations due to the high cost associated with considering the full spectrum, and because a simulation with even and odd mode numbers until n = 12, which ran for a short time, showed the same overall linear and non-linear dynamics. After a brief period of stability ( ${\sim}0.2~\mathrm{ms}$ ) the steepening pedestal that is not unstable to ideal PB modes begins to excite non-ideal PB instabilities with predominant ballooning features. The magnetic energies of the non-axisymmetric perturbations can be observed in figure 4(a). The time frame between the first vertical black line ( ${t = 0.3~\mathrm{ms}}$ ) and the vertical purple line ( ${t = 0.5~\mathrm{ms}}$ ) denotes the linear growth phase, where only the linearly unstable modes ( ${n = 6,~8,~10,~\mathrm{and}~12}$ ) grow. The time frame between the purple line and the second vertical black line ( ${t = 0.7~\mathrm{ms}}$ ) denotes the early non-linear growth phase where non-linear mode coupling excites linearly stable modes to grow (n = 2 and 4)¹¹ . The structure of the PB modes of figure 4(a) during the linear phase (at ${t = 0.4~\mathrm{ms}}$ ) is shown in figure 5, with the perturbations of density, temperature, and poloidal magnetic flux. The flux surfaces at ${\psi_\mathrm{N} = 0.95,1.00,1.05,1.10}$ are also shown as thin black lines. It is observed that the PB modes are localised very close to the separatrix, and they rotate in the electron diamagnetic direction (counter-clockwise in figure 5). Experimentally, the QCM that underlies the cross-field transport in the enhanced D-alpha H-mode is observed to travel in the electron diamagnetic direction in the lab frame [16].

**Figure 4.** Magnetic energies of the n ≠ 0 modes in a logarithmic scale in the first ${0.8~\mathrm{ms}}$ . (a) shows the simulation with ${n = 0,2,4,\dots,12}$ and (b) shows the simulation with ${n = 0,2,4,\dots,20}$ . This comparison shows that the dominant mode numbers in the linear phase are n = 10 and 12 in both cases; into the non-linear phase (not shown), the n > 12 modes are sub-dominant.
Download figure:
Standard image High-resolution image

**Figure 4.** Magnetic energies of the n ≠ 0 modes in a logarithmic scale in the first ${0.8~\mathrm{ms}}$ . (a) shows the simulation with ${n = 0,2,4,\dots,12}$ and (b) shows the simulation with ${n = 0,2,4,\dots,20}$ . This comparison shows that the dominant mode numbers in the linear phase are n = 10 and 12 in both cases; into the non-linear phase (not shown), the n > 12 modes are sub-dominant.
Download figure:
Standard image High-resolution image

**Figure 5.** Density, temperature, and poloidal magnetic flux perturbations (n ≠ 0 components) from the simulation with ${n = 0,2,4,\dots,12}$ at ${0.4~\mathrm{ms}}$ . The peeling-ballooning structure with dominant ballooning characteristics, i.e. more localised to the low field side (LFS), can be observed in both plots. Flux surfaces at ${\psi_\mathrm{N} = 0.95,1.00,1.05,}$ and 1.10 are shown as grey lines.
Download figure:
Standard image High-resolution image

**Figure 5.** Density, temperature, and poloidal magnetic flux perturbations (n ≠ 0 components) from the simulation with ${n = 0,2,4,\dots,12}$ at ${0.4~\mathrm{ms}}$ . The peeling-ballooning structure with dominant ballooning characteristics, i.e. more localised to the low field side (LFS), can be observed in both plots. Flux surfaces at ${\psi_\mathrm{N} = 0.95,1.00,1.05,}$ and 1.10 are shown as grey lines.
Download figure:
Standard image High-resolution image

The initial growth phase is started by an n = 12 mode in these simulations and closely followed by the growth of the n = 10 mode. The growth rates of these modes are very similar, roughly ${\gamma_{n = 10,12}\approx5\times10^4\ \mathrm{s}^{-1}}$ . A separate simulation including all even toroidal mode numbers until n = 20 has been produced, and the magnetic energies of the n ≠ 0 modes are shown in figure 4(b). In this way, we confirm that the fastest growing mode is the n = 12, both when the highest toroidal mode number in the system is ${n_\mathrm{max} = 12}$ and when it is ${n_\mathrm{max} = 20}$ . The toroidal harmonics with n > 14 grow at a slower rate than the n = 12 because diamagnetic stabilisation is stronger for modes with a larger toroidal mode number. Additionally, they remain sub-dominant well into the non-linear phase (not shown). Quadratic non-linear mode coupling takes place during the early non-linear phase ( ${t = [0.5,0.7]~\mathrm{ms}}$ ). The non-linear mode coupling that gives rise to the excitation of linearly stable modes has been described in JOREK simulations in [46].

During the linear growth phase (at ${t = 0.5~\mathrm{ms}}$ ), the poloidal velocity of the linearly unstable PB modes ( ${n = 6,8,10,12}$ ) in the outboard midplane is ${v_\mathrm{mode,pol}\approx11~\mathrm{km}\,\mathrm{s}^{-1}}$ (the positive sign indicating electron diamagnetic direction). Resistive and ideal ballooning modes have been identified as travelling in the laboratory frame with the following velocities,

$\begin{align} \,\mathrm{Resistive\!\!:}& \quad v_{\mathrm{mode},pol} = \,\,v_\mathrm{ExB} + v_{\parallel,\mathrm{pol}} , \end{align} \tag{ 1 }$

$\begin{align} \qquad\qquad\qquad \mathrm{Ideal\!\!:}& \quad v_{\mathrm{mode},pol} = \,\,v_\mathrm{ExB} + v_{\parallel,\mathrm{pol}} + v^*_\mathrm{i}/2, \end{align} \tag{ 2 }$

where ${v^*_\mathrm{i} = \nabla p_\mathrm{i}/(e n_e B)}$ is the ion diamagnetic velocity, and ${v_{\parallel,\mathrm{pol}}}$ is the poloidal projection of the parallel velocity [47]. At ${t = 0.5~\mathrm{ms}}$ , equation (1) results in ${v_\mathrm{mode,pol}\approx10~\mathrm{km}\,\mathrm{s}^{-1}}$ , and equation (2) results in a poloidal velocity of roughly ${v_\mathrm{mode,pol}\approx4~\mathrm{km}\,\mathrm{s}^{-1}}$ . Namely, ${v_\mathrm{mode,pol}}$ for these modes is closer to that of resistive ballooning modes than to ideal ballooning modes. The following subsection is devoted to directly studying the influence of plasma resistivity on the growth rates of the non-axisymmetric perturbations.

4.2. Rigidly scanning the resistivity

In this subsection, we freeze the axisymmetric profiles at ${t = 0.5~\mathrm{ms}}$ and change the resistivity to understand its influence on the stability of the non-ideal PB modes of all even mode numbers up to n = 12. For this scan we consider several multiplication factors of the nominal resistivity (nominal resistivity on axis is ${\eta_\mathrm{axis} = 6.6\times10^{-8}~{\Omega \textrm{m}}}$ , and it follows the Spitzer temperature dependency), ${0.5,~1.5,~3.0,~6.0,~12.0,~24.0,~48.0}$ . Varying the resistivity results in a change in the growth rate of the linearly unstable modes, as shown in figure 6. The x-axis represents the resistivity (in ${{\Omega \textrm{m}}}$ ) and the y-axis corresponds to the growth rate (in $\mathrm{1/s}$ ); both axes are plotted in logarithmic scales. Different colours and symbols represent different toroidal mode numbers.

**Figure 6.** The growth rates of different toroidal mode numbers in the resistivity scan. Increasing the resistivity leads to the destabilisation of resistive PB modes. At large resistivity, the growth rates follow a ${\eta^{1/3}}$ power law.
Download figure:
Standard image High-resolution image

**Figure 6.** The growth rates of different toroidal mode numbers in the resistivity scan. Increasing the resistivity leads to the destabilisation of resistive PB modes. At large resistivity, the growth rates follow a ${\eta^{1/3}}$ power law.
Download figure:
Standard image High-resolution image

The response of the PB modes to the changes in resistivity displays several noteworthy characteristics. At nominal resistivity, ${n = 2~\mathrm{and}~4}$ are linearly stable and become linearly unstable at 3.0 and 1.5 times the nominal resistivity, respectively (i.e., the said modes no longer require non-linear mode coupling to grow). Conversely, at half of the nominal resistivity, the n = 6 mode becomes linearly stable. This is not surprising given the fact that we know (from the MISHKA ideal MHD simulations) that ideal PB modes are not unstable for the considered profiles and, as such, at arbitrarily low resistivity, all toroidal mode numbers ought to be stable. At high resistivity, the resulting scaling is that of resistive interchange modes, as proposed by Furth et al (which breaks down when η is too large) [48]. The fact that we only recover the scaling at high η is expected, as we are simulating instabilities that cross a threshold value of resistivity to become destabilised and, therefore, below and near the threshold, it is not only the resistivity that influences the growth rates. Recent work with the M3D-C1 code has also studied how finite resistivity affects peeling-ballooning modes, but in the context of spherical tokamaks [49].

4.3. Importance of extended MHD

Simulations without diamagnetic effects have been performed in order to understand their influence on the resistive instabilities described in the previous section. This is done for simulations with only one toroidal harmonic present, n = 8, and for different applied heating powers. The nominal heating power considered (denoted $P_\mathrm{heat,0}$ ) corresponds to the simulation of section 4.1, and it is multiplied by factors of 1.0, 1.25, 1.63, and 2.51 to understand the effect of a temperature pedestal that grows at different rates. It is observed that the simulations that include $v^*_i$ have fundamentally different non-axisymmetric dynamics with respect to the simulations that neglect the diamagnetic effects. Figure 7 shows the evolution of the n = 8 magnetic energy in the logarithmic scale of simulations (a) with and (b) without diamagnetic effects, at four different multiplication factors of ${P_\mathrm{heat,0}}$ .

**Figure 7.** The evolution of ${E_\mathrm{mag,n = 8}}$ for simulations (a) with and (b) without ion pressure gradient-driven diamagnetic flows. Four different input heating powers are considered. Increasing heating power in simulations that include diamagnetic flows shows an important stabilisation of the n = 8 PB mode. When neglecting the diamagnetic effects, increasing heating power causes the unstable PB mode to grow even faster due to the steeper pressure profiles.
Download figure:
Standard image High-resolution image

**Figure 7.** The evolution of ${E_\mathrm{mag,n = 8}}$ for simulations (a) with and (b) without ion pressure gradient-driven diamagnetic flows. Four different input heating powers are considered. Increasing heating power in simulations that include diamagnetic flows shows an important stabilisation of the n = 8 PB mode. When neglecting the diamagnetic effects, increasing heating power causes the unstable PB mode to grow even faster due to the steeper pressure profiles.
Download figure:
Standard image High-resolution image

Increasing heating power causes a steepening of the edge temperature profile and, therefore, of the pressure at the plasma edge. For the simulations that include diamagnetic effects, the edge ${E_r}$ well at the pedestal becomes deeper with steeper pressure profiles. In these simulations, the PB modes become stabilised by the diamagnetic drift together with the E_r (and its shear) [33, 50, 51]. For the simulations without diamagnetic effects, the pedestal steepens, but E_r does not change. In figures 7(a) and (b), the heating power of the different simulations is changed at ${t = 0.33~\mathrm{ms}}$ to the values shown in the legend. The simulations that consider diamagnetic effects observe ${\gamma_{n = 8}}$ to decrease as the heating power is increased. On the other hand, the simulations that neglect diamagnetic effects result in an increase of ${\gamma_{n = 8}}$ with increasing heating power. The results presented in this section emphasise the importance of including the diamagnetic flows for simulations of PB modes.

As mentioned at the beginning of this section, the pedestal profiles for the nominal heating power at these stages is stable for ideal peeling-ballooning modes. Higher pedestal pressure and/or edge current densities are required in order to reach a type-I ELM unstable scenario. The access to a type-I ELM unstable scenario appears to be closed without the inclusion of the two-fluid diamagnetic effects. Indeed, this result was previously reported in [52] in the context of obtaining repetitive ELM cycle simulations, and was extended as a requirement to simulate type-I ELM cycles in [22]. The following subsections are devoted to describing the fully non-linear phase (during which the n ≠ 0 modes interact with each other and with the n = 0 axisymmetric background) of the simulation with ${n = 0,2,4,\dots,12}$ at nominal heating and resistivity in the presence of diamagnetic effects, i.e. the continuation of figure 4(a).

4.4. Non-linear phase

For the simulation that includes ${n = 0,2,4,\dots,12}$ , i.e. figure 4(a), the magnetic and kinetic energies of the non-axisymmetric modes are shown in figure 8 in linear scale for ${10~\mathrm{ms}}$ of simulation time (a) and (c) and in logarithmic scale for the first ${2~\mathrm{ms}}$ (b) and (d).

**Figure 8.** Magnetic and kinetic energies of the non-axisymmetric modes in linear scale for ${10~\mathrm{ms}}$ of simulation time (top) and in logarithmic scale for the first ${2~\mathrm{ms}}$ of simulation time, during which the linear growth phase takes place.
Download figure:
Standard image High-resolution image

The linear growth phase gives way to the early non-linear growth phase until the amplitude of the perturbations becomes large with respect to the background plasma, at which point the non-linear phase begins. During the latter, a dynamic interplay between the n ≠ 0 modes and the background plasma determines the instantaneous profiles observed in the simulations. Due to the persistent PB modes and the lack of clear cyclical dynamics, the dynamics observed can be characterised as peeling-ballooning turbulence.

4.5. Filamentary transport

The non-axisymmetric time evolution of the ${\varphi = 0}$ outer midplane pressure gradient and the inner/outer divertor incident power are shown in figures 9(a) and (b), respectively, for ${10~\mathrm{ms}}$ of simulation time. The incident power is defined as

$\begin{align*} P_\mathrm{div} = \int_0^{2\pi}\int_{\ell_0}^{\ell_\mathrm{max}}q(t,\ell,\varphi) R\,d\ell\,d\varphi, \end{align*}$

where ${q(t,\ell,\varphi)}$ is the heat flux at a given time in the divertor location $\ell$ at the toroidal angle ϕ, and R is the major radius. The colour map indicates that the pressure gradient and, therefore, the pressure profile in the outermost edge of the plasma is rapidly fluctuating. The corresponding fluctuations are governed by the non-axisymmetric modes that regulate the pedestal to fluctuate about a mean value, i.e. the PB turbulence. The incident power onto the divertors does not have characteristic spikes, but rather displays a quasi-continuous heat deposition, which is qualitatively similar to the QCE regime or the EDA H-mode in AUG. These observations allow us to distinguish our simulations as different from type-III ELMs, which appear as distinct events.

**Figure 9.** The time evolving outer midplane edge pressure gradient at ${\varphi = 0}$ in colour scale (a), and inner/outer divertor incident power (b). The varying pressure profile is caused by quasi-continuous outward transport created by non-ideal peeling-ballooning modes. The position of the separatrix is shown as a white line.
Download figure:
Standard image High-resolution image

Figure 10 shows the ${\varphi = 0}$ outer midplane pressure in a colour map with logarithmic scale for a reduced time window of ${0.4~\mathrm{ms}}$ , which is chosen between 4.8 and ${5.2~\mathrm{ms}}$ , in order to show the dynamics of filamentary structures originating from slightly inside the separatrix (roughly ${2~\mathrm{cm}}$ ) and travelling outwards. The y-axis is the major radius, and the separatrix position is represented with a white line. The plasma blobs that travel outwards result from the resistive PB modes that are aligned to the magnetic fields and are moving in the electron diamagnetic direction.

Resistive peeling-ballooning modes that are destabilised below the ideal PB stability boundary regulate the pressure gradient about ${-250~\mathrm{kPa}\,\mathrm{m}^{-1}}$ . This is made clearer with pressure and pressure gradient profiles taken in the representative time frame of ${t = [4.0,6.0]~\mathrm{ms}}$ together with time-averaged profiles in black shown in figures 11(a) and (b). The fluctuating profile in the last ${\sim}7\%$ of the confined region is clearly visible in figure 11(a). Unlike the axisymmetric simulation (figure 3), the profiles do not show a monotonic increase in pedestal top. The oscillatory nature of figures 11(a) and (b) is caused by the quasi-continuous transport incited by the resistive PB modes, and these modes are the reason why the pedestal top cannot grow further. This is why these simulations feature only small ELMs and not a mixed regime with small ELMs and type-I ELMs.

**Figure 11.** Pressure (a) and pressure gradient (b) profiles in the time window ${t = [4.0,6.0]~\mathrm{ms}}$ together with a time-averaged profile in black. The time-averaged profile shows a 'staircase' structure with a large pressure gradient in the vicinity of the separatrix.
Download figure:
Standard image High-resolution image

**Figure 11.** Pressure (a) and pressure gradient (b) profiles in the time window ${t = [4.0,6.0]~\mathrm{ms}}$ together with a time-averaged profile in black. The time-averaged profile shows a 'staircase' structure with a large pressure gradient in the vicinity of the separatrix.
Download figure:
Standard image High-resolution image

Taking the time-varying temperature fluctuations at a single point in the steep gradient region in the outer midplane, ${(R,Z) = (2.14,0.06)}$ , a spectrogram is performed. As a result, a dominant frequency in the range of $20{-}40~\mathrm{kHz}$ is found, as can be observed in figure 12. A type-II ELMy H-mode in AUG with high triangularity and close to double null reported in [7] was described as having an electron pressure gradient oscillating about ${\sim}150~\mathrm{kPa}\,\mathrm{m}^{-1}$ . The oscillating ${\nabla p_e}$ was reportedly caused by MHD modes, which were associated with electromagnetic fluctuations observed in a wide radial extent peaking at a frequency range of $30{-}50~\mathrm{kHz}$ . Both observations hint at qualitative similarities to the simulation results described in this section. Nevertheless, it must be noted that the present simulations were performed in a different magnetic configuration, i.e., with low triangularity and far from double null. Therefore, dedicated comparisons need to be performed in the future to produce quantitative comparisons between experiments and simulations. In particular, such comparisons will have to include variations of the plasma shape.

**Figure 12.** Time evolving (a) and averaged (b) frequency spectrogram of the temperature fluctuations in ${(R,Z) = (2.14,0.06)}$ . Dominant frequencies in the range of $20{-}40~\mathrm{kHz}$ can be observed in both cases.
Download figure:
Standard image High-resolution image

4.6. Divertor heat deposition

To show the quasi-continuous exhaust caused by the non-ideal peeling-ballooning modes excited near the separatrix, the electron temperature at the inner and outer divertor targets is plotted in figures 13(a) and (b), respectively. The target electron temperature is considered to be half of the plasma temperature, and it is plotted for $10~\mathrm{ms}$ of simulation time. The inner divertor target has a lower target temperature than the outer divertor. Similarly, the incident power to the inner divertor is lower than to the outer divertor, as seen in figure 9(b); this is consistent with calorimetry measurements during AUG discharges in the QCE regime [53]. There is a slight increase in the maximum target temperature (particularly visible in the outer target) as time progresses. This is due to the chosen heat source in the confined region, which slowly increases the thermal energy content inside the separatrix. At any given time point, the heat deposition does not show significant variations in the toroidal direction. In other words, the heat deposition is roughly axisymmetric. It is important to note that the present simulations use only a simplified SOL transport model and, as such, the obtained heat distribution between targets will not necessarily reflect experimental observations. A more advanced SOL/divertor model is being developed in JOREK [45] and it will be used for future simulations.

**Figure 13.** Time evolution of (a) the inner and (b) outer target electron temperature ( ${T_e = T/2}$ ) caused by the resistive PB modes excited at the very edge of the plasma. The inner target has a lower temperature than the outer target as well as a lower incident power.
Download figure:
Standard image High-resolution image

**Figure 13.** Time evolution of (a) the inner and (b) outer target electron temperature ( ${T_e = T/2}$ ) caused by the resistive PB modes excited at the very edge of the plasma. The inner target has a lower temperature than the outer target as well as a lower incident power.
Download figure:
Standard image High-resolution image

5. Two simple paths to type-I ELMs

The simulations presented in the previous section showed that the small ELM regime under consideration is sustained by quasi-continuous transport caused by resistive peeling-ballooning modes. Based on such simulations, two different ways to transition from the small-ELM dominant scenario towards a type-I ELMy H-mode are investigated. The first approach is to increase the heating power, and it is described in section 5.1. The second approach results from decreasing the separatrix density while maintaining the pedestal top density, and it is detailed in section 5.2. Both paths lead to a transition to a type-I ELMy H-mode resulting from a reduction in the drive of the resistive peeling-ballooning modes that underlie the small-ELM regime (reduction of the resistivity when ${P_\mathrm{heat}}$ is increased, and reduction of the pressure gradient near the separatrix when ${n_\mathrm{e,sep}}$ is decreased), together with an increase of stabilising effects (increase of diamagnetic and plasma flows in both cases and, additionally, of the edge current density when ${n_\mathrm{e,sep}}$ is reduced).

5.1. Increasing heating power

The nominal heating power, ${P_\mathrm{heat,0} = 6.2\times10^{-6}}$ in JOREK units, corresponds to the simulation shown in figure 8. The magnetic energies of the non-axisymmetric perturbations for the first ${7~\mathrm{ms}}$ of the simulation are shown in figure 14(a). In the subsequent sub-figures, the heating power is progressively increased¹² in small steps (the excess heating power is included from the beginning of each simulation). In figure 14(b), the non-axisymmetric activity remains strong enough to prevent the pedestal top from growing (not shown), and the n = 8 mode appears to be more dominant than in the simulation with nominal heating power. With the next increase of $P_\mathrm{heat}$ , figure 14(c), the dominance of the n = 8 mode becomes even stronger, and there is a transient phase where it hosts most of the total non-axisymmetric magnetic energy, ${\Sigma_{n\gt0}} E_{\mathrm{mag}, n}$ . In this case, the non-ideal PB modes become completely stabilised after roughly $10~\mathrm{ms}$ (not shown). Finally, in figures 14(d) and (e), the dominance of the n = 8 mode is almost complete, and ${\Sigma_{n\gt0}} E_{\mathrm{mag}, n}$ is reduced until complete stabilisation in progressively shorter time scales.

**Figure 14.** Magnetic energies of the non-axisymmetric perturbations for five different values of input heating power. The applied heating power increases progressively from (a) to (e). As a result of the increasing heating power, the magnetic energies of the n > 0 perturbations become completely suppressed if sufficient additional heating power is considered.
Download figure:
Standard image High-resolution image

The stabilisation of the resistive PB modes at higher heating power can be attributed to the fact that the pedestal top temperature increases which, in turn, decreases the local resistivity ( ${\eta\propto T^{\,-3/2}}$ ). The latter results in a decrease in the linear growth rates of the resistive PB modes, as shown in section 4.2. Less linearly unstable modes may then be stabilised through the action of diamagnetic and plasma flows more easily. This mode stabilisation allows the pedestal to build-up and eventually gives rise to a type-I ELM crash (the physical effects involved in such bifurcation are discussed and shown in more detail later in this section).

To further understand what governs the transition from small ELMs to type-I ELMs, the radial electric field at the outboard midplane is averaged between $1.0{-}2.0~\mathrm{ms}$ and it is shown in figures 15(a)–(e). It is interesting to note that the three scenarios where the small ELMs become stabilised, (c)–(e), have deeper radial electric field wells than the two cases that sustain the small ELMs, (a) and (b). It is worth pointing out that even for the lowest heating power, the instantaneous radial electric field profiles at the outer midplane are often deeper than ${E_r\sim-15~\mathrm{kV}\,\mathrm{m}^{-1}}$ , a representative value that has been associated to the L-H transition in AUG [54, 55]. This can be seen in figure 15(a), which shows in grey the instantaneous profiles used to obtain the time-averaged profile (black). The scans with and without diamagnetic effects shown in section 4.3, together with the observations presented so far in this section, indicate that small ELMs are susceptible to a stabilising effect from larger diamagnetic drifts and the resulting deeper radial electric field well.

**Figure 15.** Time-averaged profiles of the outboard midplane radial electric field for five different values of input heating power. The applied heating power increases progressively from (a) to (e).
Download figure:
Standard image High-resolution image

In order to directly show the bifurcation from small ELMs to type-I ELMs, the heating power of the simulation described in section 4.4 is increased at ${5.6~\mathrm{ms}}$ of simulation time up to ${P_\mathrm{heat} = 6.6\times 10^{-6}}$ in JOREK units, i.e., the heating power corresponding to figure 14(e). As a result of the heating power increase, the resistive PB modes start to become smaller in amplitude and their radial extent starts to reduce. This process takes roughly ${4~\mathrm{ms}}$ to complete, and as it takes place, a steeper pedestal is allowed to form. Ultimately, the steepening pedestal crosses the type-I ELM stability boundary and an ELM with dominant toroidal mode numbers n = 2 and 4 is excited. The process described here can be evidenced in figures 16(a) and (b), which respectively show the ${\varphi = 0}$ outboard midplane pressure gradient and the power incident on the inner and outer divertors. The small jump at ${t\approx12~\mathrm{ms}}$ takes place because the parallel heat conductivity has been increased roughly to the Spitzer–Härm values in this simulation, but the onset of the type-I ELM takes place regardless of the said change. The incident power that reaches the divertor targets is significantly increased when the type-I ELM crash appears. It is comparatively clear that the small ELMs cause much weaker heat fluxes to the divertor targets. To directly show the influence of the increased heating power onto the resistive PB modes that cause small ELMs, the plasma pressure in real space, together with the position of the separatrix, are plotted for a restricted time frame between t = 5 and $10~\mathrm{ms}$ in figure 17. The expelled filaments after the heating power increase seem to have smaller amplitudes, and they travel for shorter distances. They eventually disappear completely.

**Figure 16.** The time evolution of the outboard midplane pressure gradient at ${\varphi = 0}$ (a) upon increasing the heating power at ${5.6~\mathrm{ms}}$ . The power incident on the inner and outer divertor targets (b) resulting from small ELMs ( $t\lt10~\mathrm{ms}$ ) and from a type-I ELM ( $t\gtrsim19~\mathrm{ms}$ ).
Download figure:
Standard image High-resolution image

**Figure 17.** The time evolution of the outboard midplane pressure at ${\varphi = 0}$ (a) upon increasing the heating power at ${5.6~\mathrm{ms}}$ . The filaments expelled from the confined region become weaker when the heating power is increased. Ultimately, they completely disappear and the pedestal is able to grow further.
Download figure:
Standard image High-resolution image

In section 2, we discussed the response of AUG small ELMs at high separatrix density, low triangularity and high edge safety factor towards increasing heating power. Namely, in such experiments the small ELMs can become suppressed with sufficiently high additional heating power. The situation is similar for the EDA H-mode. In both cases, the role of plasma shaping in the suppression of small ELMs/QCM seems to be pivotal. In our simulations, which feature high ${n_{e,\mathrm{sep}}}$ , low triangularity, and high ${q_{95}}$ , the transition from a regime dominated by small ELMs towards a type-I ELM is obtained by suddenly increasing the input heating power in small ELM simulations. The small ELMs start to weaken and the filaments formed by the small ELMs are gradually reduced in amplitude until disappearing completely. ${p_\mathrm{ped}}$ rises with increasing ${P_\mathrm{heat}}$ due to the increase of ${T_\mathrm{ped}}$ . ${n_{e,\mathrm{ped}}}$ remains unchanged in the first few milliseconds after the heating power was increased; it only starts rising when the particle transport by small ELMs becomes significantly reduced. Additionally, due to the larger ${\nabla p}$ , the radial electric field well at the plasma edge deepens, and the bootstrap current density starts to rise. Taking 1 ms time-averages, the outer midplane profiles are tracked during the pure small ELMs phase (from ${2.6~\mathrm{ms}}$ until ${5.6~\mathrm{ms}}$ ) and during the transition phase, where the resistive PB modes start to disappear (from ${5.6~\mathrm{ms}}$ to ${9.6~\mathrm{ms}}$ ) and are plotted in figure 18.

**Figure 18.** Evolution of the time-averaged ${\varphi = 0}$ outboard midplane pressure (a), radial electric field (b), and toroidal current density (c). The time averaging is done for 1 ms intervals. The three time-averaged profiles before the heating power is increased have a shallow E_r, and the four profiles after ${P_\mathrm{heat}}$ show systematically higher ${p_\mathrm{ped}}$ and a deeper E_r well (as indicated by the black arrow for the latter).
Download figure:
Standard image High-resolution image

Three time-averaged profiles correspond to the original small ELM phase and show a roughly constant $p_\mathrm{ped}$ , a weak E_r at the edge, and low toroidal current density. On the other hand, the four profiles at higher heating power show systematically higher $p_\mathrm{ped}$ , a deeper E_r well, and a broader and higher edge current density. The last time-averaged profiles in grey, with the highest $p_\mathrm{ped}$ , ${\mathrm{min}(E_r)\approx-28~\mathrm{kV}\,\mathrm{m}^{-1}}$ and a high toroidal current density, is taken during a phase that mostly has suppressed the resistive PB modes. The destabilising effect of the increased ${\nabla p}$ is tied to the stabilising effects of the deepening of the E_r well and of the toroidal current density. The additional stabilising effect of the edge resistivity decreasing as the pedestal top temperature increases, is also fundamental at this stage. Therefore, the disappearance of the resistive PB modes appears to be due to the decreasing local η together with the deepening E_r and the $-j_\varphi$ increasing.

5.2. Decreasing separatrix density

The previous section detailed the bifurcation from a small ELM-dominant regime to a type-I ELM, by means of increasing the heating power. Another path towards type-I ELMs, starting from a small ELM regime, is to decrease the separatrix density. In the experiment, this can be achieved by reducing, or completely removing, the particle source given by a gas puff (replacing the particle flux by means of cryogenic deuterium pellet injection can keep $n_{e,\mathrm{ped}}$ approximately unchanged). Indeed, it has been reported that separatrix densities below ${\sim}0.35 \times n_{GW}$ do not host dominant small ELM phases [9].

The simulations presented thus far have a separatrix density of ${n_{e,\mathrm{sep}}\approx3\times10^{19}~\mathrm{m^{-3}}}$ . To reduce the separatrix density in our simulations, we reduce the edge density source. In order to maintain ${n_{e,\mathrm{ped}}}$ unchanged, we increase the core density source. This has been achieved experimentally by reducing gas fuelling while, at the same time, increasing pellet fuelling to incite a change from a small-ELM dominant regime to a type-I ELM dominant regime [8]. As a result of the lower ${n_{e,\mathrm{sep}}}$ , the pressure gradient is locally reduced, and it is shifted slightly inwards overall; such a response of the pressure profile position is also observed in experiments [56, 57]. The decrease of the separatrix density (at fixed ${n_{e,\mathrm{ped}}}$ ) also causes a deeper ${E_r}$ well (because ${E_r\propto1/n_e}$ ) and a higher bootstrap current density (because the density gradient increases). An important caveat must be mentioned: many physical effects related to the pedestal position (particularly related to the penetration of neutrals) are not included in the JOREK model used for these simulations and, therefore, the qualitative agreement will likely not translate to a quantitative agreement at this stage.

A new simulation is then set-up with lower ${n_{e,\mathrm{sep}} = 2\times10^{19}~\mathrm{m^{-3}}}$ . The non-zero toroidal modes are included after ${0.1~\mathrm{ms}}$ (exactly the same time as the small ELM simulations described in section 4). The magnetic energies of the high and low ${n_{e,\mathrm{sep}}}$ simulations are shown in figures 19(a) and (b), respectively. The linear phases are similar between the two cases, with growing $n = 8,10,12$ PB modes. The simulation with low $n_{e,\mathrm{sep}}$ (b) deviates as the n < 10 and n = 12 become completely stabilised, and the n = 10 only reaches small amplitudes and does not affect the n = 0 background. The stabilisation of the modes with higher toroidal mode numbers leads to a single toroidal mode number with a very small amplitude that does not cause any changes to the background plasma.

**Figure 19.** Magnetic energies of the non-axisymmetric modes in the first 4 ms of simulation time for (a) small ELMs at high ${n_\mathrm{sep}({\sim}3\times10^{19}~\mathrm{m^{-3}})}$ and (b) their response to lower separatrix density ( ${\sim}2\times10^{19}~\mathrm{m^{-3}}$ ), shown in logarithmic scale .
Download figure:
Standard image High-resolution image

The simulation at lower separatrix density, figure 19(b), sees the pedestal evolve, but it does not continue until a type-I ELM is reached to save computing time. Time-averaged outer midplane profiles of the pressure gradient and radial electric field are displayed in figure 20 for the small ELMs with ${n_{e,\mathrm{sep}}\approx3\times10^{19}~\mathrm{m^{-3}}}$ (a) and (c), and for the lowered separatrix density case ( ${n_{e,\mathrm{sep}}\approx2\times10^{19}~\mathrm{m^{-3}}}$ ) (b) and (d). The profiles are averaged over $0.1~\mathrm{ms}$ during the linear growth phase, ${\lt}0.7~\mathrm{ms}$ . The lower $n_\mathrm{sep}$ causes an inward shift of $\nabla p$ and, particularly, a smaller pressure gradient near the separatrix. This additionally allows for a deeper E_r well and an increase in the bootstrap current density (not shown). The reduced drive of the $\nabla p$ near the separatrix together with the stabilising influence of the deeper E_r and the higher $-j_\varphi$ cause the small ELMs to become completely stabilised.

**Figure 20.** The time-averaged outboard midplane profiles of the pressure gradient and E_r during the linear growth phase for a nominal separatrix density (a) and (c), and for a lower separatrix density (b) and (d).
Download figure:
Standard image High-resolution image

6. Conclusions

H-mode operation without large type-I ELMs is an imperative requirement for ITER in high-performance conditions. To this end, naturally ELM-free H-modes and ELM mitigated/suppressed regimes are considered and actively investigated. Several such alternatives are under investigation in AUG. Two of them are the QCE regime and the enhanced D-alpha H-mode. Both can be operated completely without type-I ELMs. The pedestal is limited by small ELMs for the QCE regime and by a QCM for the EDA H-mode. These transport mechanisms quasi-continuously expel heat and particles from the confined region. It is presently unclear whether it will be possible to operate such regimes in ITER. Simulations performed with JOREK, which show several key features of small ELMs, are presented in this paper. Resistive peeling-ballooning modes near the separatrix are identified as the transport mechanism underlying such small ELMs.

The pedestal build-up is modelled at a fixed pedestal width with stationary diffusion coefficients and sources, following the approach described in [22]. Under appropriate conditions, resistive peeling-ballooning modes that regulate the pedestal below the ideal PB stability boundary are observed. The necessary conditions to sustain sufficient outward transport by small ELMs are primarily determined by the separatrix density and the input heating power. In particular, simulations with high $n_{e,\mathrm{sep}}$ and low heating power observe phases (longer than $10~\mathrm{ms}$ ) with quasi-continuous outward transport that prevent the pedestal from reaching a type-I ELM unstable scenario. The resistive nature of such PB modes is determined by the fact that they appear below the ideal PB stability boundary, because their growth rates are largely reduced/enhanced by decreasing/increasing resistivity, and because their poloidal mode velocities are measured to be close to that expected for resistive modes. The fluctuations caused by the resistive PB modes are measured to be in the frequency range of $20{-}40~\mathrm{kHz}$ , which is comparable to those measured experimentally for small ELMs [7, 29] and for the QCM in EDA H-mode [16, 29].

An important ingredient required, in order to properly simulate these resistive PB modes, is the inclusion of diamagnetic effects, which (in the simulations) allow the radial electric field well to develop in the pedestal region. In the absence of diamagnetic effects, it is not possible to stabilise the small ELMs by increasing the heating power. In contrast, when diamagnetic effects are included, the small ELMs become completely stabilised and the plasma state moves to a type-I ELMy H-mode by increasing ${P_\mathrm{heat}}$ . Similarly, decreasing the separatrix density completely stabilises the small ELMs if diamagnetic effects are included. Another important effect that should be included when modelling these instabilities is the bootstrap current density, because it has a stabilising influence on high-n peeling-ballooning modes. At the moment, JOREK evolves the bootstrap current density through the Sauter formula [34, 35], as explained in section 3. However, Sauter's expression is known to be inaccurate depending on the parameter regime, particularly at high collisionality [58]. Therefore, an improvement of the bootstrap current density source in JOREK will be pursued in the future. Finally, the simplified resistivity with only Spitzer temperature dependency used in JOREK should be improved to include the influence of neoclassical effects and an effective main ion charge greater than unity.

Acknowledgments

This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No 101052200 — EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them. In particular, contributions by EUROfusion work packages Enabling Research (EnR), Medium Size Tokamaks (MST), Working Package Tokamak Exploitation, and TSVV Project 8 are acknowledged. The simulations presented in this work were performed on the Marconi-Fusion supercomputer operated by CINECA in Italy.

Data availability statement

The data that support the findings of this study are available upon reasonable request from the authors.

MHD simulations of small ELMs at low triangularity in ASDEX Upgrade

Article metrics

Submit

Author e-mails

Author affiliations

Author notes

ORCID iDs

Dates

Peer review information

Abstract

1. Introduction

2. Small ELMs and EDA H-mode at ASDEX Upgrade

3. JOREK simulation set-up and axisymmetric build-up

3.1. Numerical parameters and simulation set-up

3.2. Axisymmetric pedestal build-up

3.3. Limitations of the present approach

4. Non-axisymmetric simulations

4.1. Linear and early non-linear phases

4.2. Rigidly scanning the resistivity

4.3. Importance of extended MHD

4.4. Non-linear phase

4.5. Filamentary transport

4.6. Divertor heat deposition

5. Two simple paths to type-I ELMs

5.1. Increasing heating power

5.2. Decreasing separatrix density

6. Conclusions

Acknowledgments

Data availability statement

Footnotes

MHD simulations of small ELMs at low triangularity in ASDEX Upgrade

Article metrics

Submit

Share this article

Author e-mails

Author affiliations

Author notes

ORCID iDs

Dates

Peer review information

Abstract

1. Introduction

2. Small ELMs and EDA H-mode at ASDEX Upgrade

3. JOREK simulation set-up and axisymmetric build-up

3.1. Numerical parameters and simulation set-up

3.2. Axisymmetric pedestal build-up

3.3. Limitations of the present approach

4. Non-axisymmetric simulations

4.1. Linear and early non-linear phases

4.2. Rigidly scanning the resistivity

4.3. Importance of extended MHD

4.4. Non-linear phase

4.5. Filamentary transport

4.6. Divertor heat deposition

5. Two simple paths to type-I ELMs

5.1. Increasing heating power

5.2. Decreasing separatrix density

6. Conclusions

Acknowledgments

Data availability statement

Footnotes