# Nonlinear simulations of the peeling-ballooning instability of super H-modes in the HL-3 tokamak

Yiren Zhu, Wulyu Zhong, Jiaxian Li, Miao Xue, Zhengji Li, Lei Xue, Guoliang Xiao, Ruihai Tong, Yunpeng Zou, Jiao Yin, et al.

► To cite this version:

Yiren Zhu, Wulyu Zhong, Jiaxian Li, Miao Xue, Zhengji Li, et al. Nonlinear simulations of the peeling-ballooning instability of super H-modes in the HL-3 tokamak. Nuclear Fusion, 2024, 64 (9), pp.096019. 10.1088/1741-4326/ad5e94. cea-04788881

HAL Id: cea-04788881

https://cea.hal.science/cea-04788881v1

Submitted on 18 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Nonlinear simulations of the peeling-ballooning instability of super H-modes in the HL-3 tokamak

Yiren Zhu<fnref n="1" />, Wulyu Zhong<fnref n="1" />,*, Xiaolan Zou<fnref n="2" />, Jiaxian Li<fnref n="1" />, Miao Xue<fnref n="1" />, Zhengji Li<fnref n="1" />, Lei Xue<fnref n="1" />, Guoliang Xiao<fnref n="1" />, Ruihai Tong<fnref n="1" />, Yunpeng Zou<fnref n="1" />, Jiao Yin<fnref n="1" />, Anshu Liang<fnref n="1" />, Linzi Liu<fnref n="1" />, Boda Yuan<fnref n="1" />, Min Xu<fnref n="1" /> and HL-3 Team<fnref n="1" />

E-mail: zhongwl@swip.ac.cn

Received 21 November 2022, revised 25 May 2024

Accepted for publication 3 July 2024

Published 30 July 2024

<figure><img src="image_2.png" /></figure>

CrossMark

## Abstract

As the newly built tokamak in China, HL-3 will explore high-performance operation scenarios, such as super H-mode. The energy confinement and core parameters in the super H-mode can be much larger than that in the normal H-mode. Based on the pedestal simulation code EPED, the operation space of the super H-mode is obtained in HL-3. Magnetic shear decreases with increasing triangularity; consequently, a super H-mode can be achieved. The threshold of triangularity for accessing a super H-mode in HL-3 is around 0.4. By using BOUT++, a nonlinear simulation study of the pedestal instabilities in the super H-mode equilibrium is executed for the first time. As expected, the low *n* peeling mode, which can cause much of the energy loss (17%) from the pedestal region, is dominant in the super H-mode. Such a large collapse in the pedestal region would lead to a transition from super H-mode to H-mode. It is crucial to expand the parameter space of the super H-mode or mitigate the edge-localized mode (ELM) size for sustaining the super H-mode operations. The *E* × *B* velocity shear is found to play an important role in controlling the ELMs in HL-3. The small *E* × *B* velocity shear leads to a large growth rate but results in a small ELM size around the peeling boundary. The ELM size is closely related to both the growth rate of peeling-ballooning mode and the duration time of the linear phase. In contrast, a large *E* × *B* velocity shear can stabilize the instabilities near the ballooning boundary. Next, the parameter space of the super H-mode can be enlarged.

Keywords: super H-mode, ELM, EPED, *E* × *B* velocity shear, HL-3

## 1. Introduction

H-mode is the proposed operation scenario in future tokamak devices, such as international thermonuclear experimental reactor (ITER) [1] and China Fusion Engineering

Test Reactor (CFETR) [2]. The edge transport barrier formed in the H-mode is usually referred as the pedestal pressure at the pedestal top in a tokamak normally determines the global confinement and fusion performance [3]. While developing a high-performance H-mode, a new ‘super H-mode’ regime was first predicted by the EPED model and then validated by experiments on DIII-D [4]. The super H-mode experiments have led to the highest pressure and stored energy on DIII-D [5–7]. The Alcator C-Mod also broke the record of pedestal pressure within the super H-mode operation regime [5]. Moreover, the super H-mode

<figure><img src="image_3.png" /></figure>

<footnote n="1">¹ Southwestern Institute of Physics, PO Box 432, Chengdu 610041, China</footnote>

<footnote n="2">² CEA, IRFM, F-13108 Saint-Paul-lez-Durance, France</footnote>

<footnote>* Author to whom any correspondence should be addressed.</footnote>

<footnote>Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.</footnote>
regime is also beneficial to realize high-performance plasma in ITER and power plant configurations, such as ARIES-ACT1 [8].

From the super H-mode experiments and simulations on DIII-D and Alcator C-Mod, it has been found that there are two requirements for accessing the super H-mode regime: strong shaping and large $E \times B$ velocity shear [7]. As shown in [4], only for strong shaping, the pedestal pressure bifurcates into three solutions when the density increases above a critical value. The upper two bifurcations define the operation parameter space of the super H-mode. In the high pressure region, the peeling mode dominates in the super H-mode, whereas the ballooning mode dominates in the normal H-mode due to the large collisionality and the low edge current density. Very high thermal energy confinement is achieved in DIII-D with large plasma rotation. It has been illustrated that a large $E \times B$ velocity shear can cause a shift of the turbulence destabilization critical gradient [9, 10] and lead to a reduction of turbulence transport, thus improving the energy confinement.

Edge-localized modes (ELMs) repetitively occur at the plasma edge and eject large transient heat loads on the plasma facing components in the H-mode regime [11, 12]. The ELM with large amplitude can cause the transition from the super H-mode to the normal H-mode with the dominant ballooning mode due to the narrow parameter space of the super H-mode [4]. To sustain the super H-mode operation and reduce the heat load on the divertor target, it is essential to mitigate the ELMs or expand the parameter space of the super H-mode. Meanwhile, ELM mitigation will also reduce the heat loads at the divertor target [13, 14].

As a newly built tokamak in China, HL-3 is designed to explore the operation scenario with high performance to support the critical issues of ITER and CFETR [15, 16]. The super H-mode regime is desired for HL-3. However, there are many concerns in the super H-mode studies of HL-3 such as large ELMs. It is essential to explore the parameter space and ELM characteristics in the super H-mode regime of HL-3. In this paper, the super H-mode operation on HL-3 is predicted by the EPED model. The operation spaces with different densities and other key parameters, such as triangularity, are determined. We are optimistic that normal H-mode operation starts at the peeling boundary on HL-3. By combining the pedestal profiles provided by the EPED and the integrated workflow on HL-3, the equilibria of the super H-mode and the normal H-mode are generated self-consistently. A nonlinear analysis of edge stability is then conducted by the BOUT++ code. Based on the nonlinear simulation results with a scan of $E \times B$ velocity shears, the physics mechanisms on how to affect the ELM activity in the super H-mode regime are illustrated. A method to expand the parameter space of super H-modes is also presented.

The paper is organized as follows. In section 2, the pedestal code EPED is introduced, and then the super H-mode regime in HL-3 is analyzed. In section 3, the effects of $E \times B$ velocity shear and current density on ELMs in

the super H-mode regime using nonlinear simulations using BOUT++ are presented. Section 4 shows the conclusions and discussions.

## 2. Operation space of super H-modes in HL-3

A super H-mode regime is proposed as a high-performance scenario for HL-3. It is needed to perform the stability analysis by theoretical simulation studies before high plasma performance experiments are executed on HL-3. There are several codes that can predict the pedestal profiles: EPED [17], EUROPED [18], and REPED [19, 20]. In this paper, we use EPED to perform a study on the prediction of super H-modes in HL-3.

As a first-principles predictive model of the pedestal height and width, EPED was developed around 2011 [21]. The EPED model comprises three parts. The first part is an equilibrium code that can use the 0D input parameters (plasma current, toroidal magnetic field, triangularity, etc) to generate a series of equilibria. In these equilibria, the density at the pedestal top is fixed, and only the pedestal temperature is scanned. The second part is a peeling-ballooning mode (PBM) stability analysis with the scan of equilibria [22]. The third part is the ballooning critical pedestal technique, which provides the analysis of the kinetic ballooning mode (KBM). Finally, the intersection of the PBM and KBM stability boundaries is considered for the maximum values of the pedestal height and width, as shown in [21].

Achieving a super H-mode on HL-3 can contribute to the development of high-performance scenarios in ITER and other future devices. The main parameters of HL-3 are presented in table 1. Based on the results of DIII-D [23], accessing the super H-mode regime can be expected in HL-3 with a large triangularity (0.5). RMP and high torque neutral beam injection (NBI) can facilitate access to super H-mode. Both will operate in the experiments on HL-3.

EPED is performed with the parameters shown in table 1. According to the scan results of EPED, the peeling-ballooning boundary of HL-3 is shown in figure 1. The density scan results of EPED show that there is a bifurcation when $n_{ped}$ is larger than about $7 \times 10^{19} \text{ m}^{-3}$ (0.78$n_{GW}$). The upper part of this bifurcation indicates a high pedestal performance scenario, which is known as the super H-mode. In this scenario, peeling modes that are driven by edge currents are dominant, which may lead to large ELM sizes. The pedestal height increases with the pedestal density. The lower part of the bifurcation indicates a normal H-mode region with large density and high collisionality. Ballooning modes are destabilized by the dominant pressure gradient in this scenario. It can be seen that the parameter space of super H-modes is very narrow. If the pedestal height is caused to be lower than the upper bifurcation caused by the ELM crash, the scenario will transform from the super H-mode to the H-mode with high collisionality.

To explore the super H-mode operation domain, the scan of magnetic field ($B_t$) and triangularities ($\delta$) has been executed
<table><caption>**Table 1.** Main parameters of HL-3 operation scenario in this study.</caption><thead><tr><th>Main parameters</th><th>Values</th></tr></thead><tbody><tr><td>Plasma current I<sub>p</sub></td><td>1 (MA)</td></tr><tr><td>Toroidal magnetic field B<sub>t</sub></td><td>1.88 (T)</td></tr><tr><td>Major radius R<sub>0</sub></td><td>1.78 (m)</td></tr><tr><td>Minor radius a</td><td>0.6 (m)</td></tr><tr><td>Elongation κ</td><td>1.7</td></tr><tr><td>Triangularity δ</td><td>0.5</td></tr><tr><td>Pedestal density n<sub>ped</sub></td><td>2.6 (10<sup>19</sup> m<sup>-3</sup>)</td></tr><tr><td>Normalized beta β<sub>N</sub></td><td>3</td></tr><tr><td>Pedestal effective charge Z<sub>eff.ped</sub></td><td>2</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Figure 1.** PBM stability diagram in $n_{ped}$ and pedestal height of HL-3. The blue line indicates the peeling boundary, and the red line indicates the ballooning boundary. The green arrow represents the parameter space of the super H-mode operation in HL-3.</figcaption></figure>

by EPED. Figure 2(a) shows that the pedestal height is found to increase with increasing magnetic field and triangularity at fixed density. It is obvious that a small $B_t$ or $\delta$ severely limits the pedestal height. Moreover, the super H-mode regime cannot be achieved when the $\delta$ is very small (0.2). The larger $B_t$ and $\delta$ can extend the peeling boundary for a wider range of $n_{ped}$, which will also lead to a higher pedestal threshold for the bifurcation in HL-3. Indeed, whether the super H-mode can be accessed is determined by the occurrence of bifurcation at the peeling-ballooning boundary. The bifurcation represents a transition from the one solution to three solutions in EPED. Figure 3 shows a number of solutions using stability calculations with EPED as a function of triangularity. The threshold of triangularity for accessing the super H-mode in HL-3 is found to be 0.4.

For a deeper understanding of the physical differences with various triangularities, the correlations among edge magnetic shear, safety factor and triangularity are plotted in figure 4. In figure 4, pedestal densities are fixed to be $6 \times 10^{19} \text{ m}^{-3}$ for each case. It is obvious that the safety factor increases with triangularity in figure 4(b), and the magnetic shear decreases when the triangularity becomes larger (figure 4(a)). The cases of HL-3 with larger triangularities and magnetic fields are more likely to be located in the second stable region in the $s-\alpha$ parameter space due to the smaller magnetic shear. This result can be explained by the fact

that lower magnetic shear (deeper magnetic well) causes the unstable region to be narrower. Moreover, the magnetic well is proportional to the safety factor [24]. HL-3 can achieve higher pedestal height with larger density in the second stable region.

$B_t = 1.88 \text{ T}$ and $\delta = 0.5$ are chosen to perform an initial study of the super H-mode on HL-3. The Greenwald density limit of $n_{GW}$ is around $8.85 \times 10^{19} \text{ m}^{-3}$, which is derived from the current and minor radii from table 1. We selected the case of $n_{ped} = 8 \times 10^{19} \text{ m}^{-3}$ to conduct the ELM analysis in section 3.

## 3. Analysis of ELMs in the HL-3 super H-mode regime

### 3.1. Equilibria construction

The first step of performing a reliable edge stability analysis is to construct a self-consistent equilibrium [25]. Integrated simulation workflow of equilibrium construction on HL-3 follows the method developed in [26] that has been validated to be reliable in devices such as ASDEX-Upgrade [27], KSTAR [28], DIII-D [29], and EAST [30]. Figure 5 shows that the pedestal pressure and edge current in the super H-mode are larger than that in the normal H-mode. The density at the pedestal top of the super H-mode regime is $n_{ped} = 8 \times 10^{19} \text{ m}^{-3}$. ELMs occur in the super H-mode regime and will cause the collapse at the pedestal region, which will lead to the transition from super H-mode to H-mode. Moreover, the larger ELM in the super H-mode will also be harmful to the divertor target. It is essential to study the characteristics of ELMs in super H-modes. We then develop and optimize the ELM control methods in this regime to sustain the super H-mode operation.

### 3.2. Comparison of ELM activities between H-mode and super H-mode

The nonlinear analysis of different equilibria is executed by BOUT++ 3-field model [31]. Based on the peeling-ballooning model with non-ideal effects such as diamagnetic effect and resistivity, the nonlinear equations for perturbations of the vorticity $\varpi$, pressure $P$, and magnetic flux $A_\|$ are shown below:

$$ \frac{\partial \tilde{\varpi}}{\partial t} + v_E \cdot \nabla \tilde{\varpi} = B_0 \nabla_\| \tilde{J}_\| + 2b_0 \times \kappa_0 \cdot \nabla \tilde{p} + \mu_{i,\|} \partial_{\|0}^2 \tilde{\varpi} + \mu_{i,\perp} \nabla_\perp^2 \tilde{\varpi} \quad (1) $$

$$ \frac{\partial P}{\partial t} + v_E \cdot \nabla P = \chi_\| \partial_{\|0}^2 P \quad (2) $$

$$ \frac{\partial \tilde{A}_\|}{\partial t} = -\nabla_\| \Phi + \frac{\eta}{\mu_0} \nabla_\perp^2 \tilde{A}_\| - \frac{\eta_H}{\mu_0} \nabla_\perp^4 \tilde{A}_\| \quad (3) $$

$$ \tilde{\varpi} = \frac{n_0 M_i}{B_0} \left( \nabla_\perp^2 \tilde{\varphi} + \frac{1}{n_0 Z_i e} \nabla_\perp^2 \tilde{p}_i \right), \\ \Phi = \tilde{\varphi} + \Phi_0, \quad \text{\*20c} $$

$$ \text{with } P = \tilde{p} + P_0 \text{ in the above equations} \quad (4) $$
<figure><img src="image_5.png" /><figcaption>Figure 2. Pedestal height with various of (a) magnetic fields and (b) triangularities.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 3. Number of solutions using stability calculations with EPED for various triangularities. The red shadow indicates the super H-mode area.</figcaption></figure>

$$J_{\parallel} = J_{\parallel 0} - \frac{1}{\mu_0} \nabla_{\perp}^2 \tilde{A}_{\parallel}, v_E = \frac{1}{B_0} (b_0 \times \nabla_{\perp} \Phi) . \quad (5)$$

where $\nabla_{\parallel} F = B \partial_{\parallel} (F/B)$ for any $F$, $\partial_{\parallel} = \partial_{\parallel 0} + \tilde{b} \cdot \nabla$, $\tilde{b} = \nabla \tilde{A}_{\parallel} \times b_0/B$, $\partial_{\parallel 0} = b_0 \cdot \nabla$, $\kappa_0 = b_0 \cdot \nabla b_0$, ion pressure $P_i$ and total pressure $P = P_i + P_e$. The radial simulation region of the normalized poloidal flux is 0.85–1.1, and the grid resolution for accurate simulation is $n_{\psi} = 260, n_{y} = 64, n_{z} = 65$, where $y$ is the poloidal coordinate and $z$ is the toroidal angle. According to [31], the parallel thermal conductivity is set to be 0. The parallel viscosity $\mu_{i,\parallel}$ is set to be $0.1\tau_A/R^2$, the perpendicular viscosity is $\mu_{i,\perp} = 0$. The resistivity $\eta$ and hyper-resistivity $\eta_H$ are determined by the Lundquist number $S$ and hyper-Lundquist number $S_H$. Here, the resistivity $\eta = \mu_0 R v_A / S$, hyper-resistivity $\eta_H = \mu_0 R^3 v_A / S_H$, and the hyper-Lundquist number $S_H = S / \alpha_H$. In the simulations of this study, the values of $S$ and $\alpha_H$ are chosen to be $1 \times 10^8$ and $1 \times 10^{-4}$, which are relevant to modestly sized tokamaks [32]. Moreover, only $n = 0$ mode of pressure perturbation are included in the simulation. As illustrated in [33], at the second crash in the ELM crash process, the zonal pressure and zonal flow play dominant roles in determining the ELM size. This part is not well considered in our research.

For efficiency, while performing nonlinear simulations in the study of the HL-3 super H mode by BOUT++, one-fifth of the torus is primarily simulated. The results show some limitations in the simulations in terms of simulation domain. The first is that the one-fifth of the torus geometry may generate poor resolution to simulate ELMs with low toroidal number (n) peeling or ballooning modes. Moreover, nonlinear mode-mode coupling might also be limited in these simulations. As illustrated in figure 4 of [34], mode-mode coupling can exert a significant effect on the saturation of both low n and high n ELMs. The full nonlinear dynamics may not be captured by only keeping specific toroidal harmonics [35]. However, figure 6 shows that the simulation results of BOUT++ with the super H-mode equilibrium in HL-3 mentioned in section 3.1 with one-fifth, one-third, and full torus have nearly the same results and trends in the evolution of ELM size and root mean square pressure (RMSP). Although the saturation of different n modes will be influenced by various tori, the growth rate of the dominant modes and the ELM sizes may not change much. Therefore, in this study, one-fifth of the torus is simulated for efficient computer time, which has also been confirmed to be reliable in the previous studies [25, 36].

Figure 7 shows that although the duration time of the linear phase of the H-mode is longer, the amplitude of the perturbed pressure in the super H-mode is larger due to the much higher growth rate. Moreover, the linear growth rates with different toroidal numbers (n) of the H-mode and super H-mode are shown in figure 8. The dominant modes for both of these two cases are low n = 5 mode, which are consistent with the parameter space near the peeling boundary. The collapse in the pedestal region can be clearly observed in figures 9(a) and (d). The ratio of the energy loss to the pedestal stored energy in the H-mode is 8%, whereas, in the super H-mode, it is much larger (17%). Such a large ELM size is supposed to be harmful to the divertor target in future fusion devices [13, 37]. In this paper, the ELM size is calculated by $\Delta W / W_{ped}$, which represents the ratio of energy loss during an ELM crash to pedestal stored energy.

For a more detailed understanding of the ELMs in H-mode and super H-mode, the evolutions of mode structures and pedestal pressure gradients are plotted in figures 9 and 10.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Evolutions of (a) magnetic shear and (b) safety factor with different triangularities. The pedestal density is fixed as $6 \times 10^{19} \text{ m}^{-3}$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Comparisons of (a) pedestal height and (b) edge current between normal H-mode and super H-mode on HL-3.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** The comparisons of (a) ELM size and (b) RMSP with various parts of torus in the simulations of BOUT++ with the super H-mode equilibrium mentioned in section 3.1.</figcaption></figure>

Figures 9(a) and (d) show the evolution of the pedestal pressure profiles in H-mode and super H-mode. In the super H-mode, the collapse (energy loss ~17%) due to an ELM crash is much larger than that in the H-mode. Figures 9(b), (c), (e), and (f) illustrate the mode structures with different toroidal numbers of the linear and saturation phases in the H-mode and super H-mode. It is obvious in figures 9(b) and (e) that the low *n* peeling mode (*n* = 5) is dominant in both cases. After the collapse in the pedestal, the current driving peeling mode is still dominant at the pedestal region in the super H-mode,

whereas, the dominant mode in the H-mode transforms to be *n* = 10. This phenomenon can be explained by figure 10.

The size of the ELM in the super H-mode is very large, and the collapse at the pedestal region leads to the weakening of the pedestal pressure in a short time (figure 10(b)). However, in H-mode, as shown in figure 10(a), the pedestal pressure profile is not completely flattened by the ELM due to its smaller size. Some higher *n* modes are destabilized due to the radial shift of the pressure gradient after the ELM crash. Figure 10(b) clearly shows that a single ELM crash in the super H-mode can cause
<figure><img src="image_10.png" /><figcaption>Figure 7. Time evolutions of root mean square (RMS) pressure at peak pressure gradient place in the H-mode and super H-mode. $\tau_A$ represents the Alfven time.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 8. Growth rates with different toroidal number of H-mode and super H-mode. Low toroidal number mode dominates in both of them.</figcaption></figure>

much energy loss from the pedestal region. The pressure gradient is flattened, which will stabilize the higher $n$ modes even further. Figures 11 and 12 show the contour plots of perturbed pressure at the outer midplane in H-mode and super H-mode, which will help to understand the process of mode transition. In figure 11, we see that $n = 5$ mode dominates at the linear phase ($t = 300 \tau_A$), which is indicated by the one light area and the dark area in one-fifth of the torus. The light area represents the transport to the separatrix, and the dark area represents the collapse in the inner region. After the ELM crash ($t = 400 \tau_A$), the number of light areas or dark areas increases to be 2 and $n = 10$ mode dominates. It is obvious that the low $n$ mode dominates before and after the ELM crash in the super H-mode from figure 12.

In conclusion, the super H-mode is dominated by the low $n$ mode, and that the ELM caused by this instability can lead to much energy loss from the pedestal in a very short time. How to suppress or mitigate the ELM to sustain super H-mode operation will be discussed in the next section.

### 3.3. Effect of $E \times B$ velocity shear on the ELM in super H-mode

Two methods to sustain the super H-mode operation are studied in this paper. The first is mitigating the ELM size and then the transition from super H-mode to high density H-mode will not occur. The second is enlarging the parameter space of super H-modes by suppressing the ELMs. In this paper, we focus on the correlations among the ELM size, $E \times B$ velocity shear, and the type of stability boundary.

The theory of the stabilized ballooning mode (high $n$ mode) with large shear flow has been known by a series of research works [38–40]. In the study of the super H-mode in HL-3, the effect of shear flow on the peeling ballooning mode is considered in the drift-ordered fluid equations mentioned in [41].

Normally, $E \times B$ velocity shear can be modified by changing the rotation artificially [42]:

$$ \frac{d\Phi_{v0}(\psi)}{d\psi} = D_0 [1 - \tanh(D_s(x - x_0))] + C, \quad (6) $$

$D_s$ determines the amplitude of modified $E \times B$ velocity shear.

Nonlinear simulations with different $E \times B$ velocity shears are executed by BOUT++ three-field module, whereas other parameters are fixed. As mentioned in [43], the $E \times B$ velocity shear not only affects the growth rate of peeling-ballooning instability but also plays an important role in the saturation phase of the ELMs. Indeed, $E \times B$ velocity shear may has different effects on the peeling mode and ballooning mode, and two equilibria are chosen to perform this study. One is the super H-mode case mentioned in figure 5, and the other is the equilibrium near the ballooning boundary, which is shown in figure 19.

#### 3.3.1. Effect of $E \times B$ velocity shear on the ELM at peeling boundary.
First, the linear growth rates of different $n$ and the evolutions of RMSP with various $E \times B$ velocity shears are shown in figure 13. The case with small shear has a larger growth rate but smaller RMSP size. In figure 14, the case with a smaller $E \times B$ velocity shear has a lower ELM size. In addition, the instabilities in figure 14(a) saturate faster than those in figure 14(b). For a further understanding of this phenomenon, a scan with different $E \times B$ velocity shear is worked out. In Figure 15, although the growth rate of the most unstable instability is well decreased with larger $E \times B$ velocity shear, the ELM size becomes larger due to the longer time of linear phase $t_l$. With the equilibrium on the peeling boundary, if the $E \times B$ velocity shear cannot totally suppress the ELM, the ELM size may be larger. It is obvious that the ELM size cannot be determined only by the growth rate or the duration of the linear phase. Indeed, the ELM size is closely correlated with
<figure><img src="image_12.png" /><figcaption>**Figure 9.** (a) Evolution of pedestal pressure, (b) mode structure of linear phase, (c) mode structure of saturation phase for H-mode; (d) evolution of pedestal pressure, (e) mode structure of linear phase, (f) mode structure of saturation phase for super H-mode.</figcaption></figure>

the saturation level. The nonlinear interaction including $n = 0$ generation will be essential to determine the ELM size [33, 44]. Further study may be executed according to the optimized model in [44], which can simulate the full annular tokamak with the consideration of $n = 0$ mode for $\phi$ and $A_\|$ better.

For further insights into the correlation among the growth rate, ELM size and other parameters, the evolutions of instabilities with different $E \times B$ velocity shear are studied in more detail. In the case of $D_s = 10$, only the instabilities with low frequency dominate during the saturation phase, whereas instabilities with several frequencies dominate in the saturation phase for the case with $D_s = -10$ (figure 16). This

indicates that a larger $E \times B$ velocity shear can suppress varies of instabilities during the ELM crash. In this situation, the low frequency modes that drive the ELM will grow to have an extended size. However, when the $E \times B$ velocity shear is small, varying instabilities are expected to remain for the whole process of ELM crash. These instabilities with wide frequencies may limit the growth of ELM; therefore, the ELM size in the saturation phase is smaller in the case with lower $E \times B$ velocity shear. Moreover, we show that the higher $n$ mode ($n = 10$) dominates at saturation phase from figure 16(b). According to figures 15 and 17, we can conclude that the $E \times B$ velocity shear has many effects on the growth rate and
<figure><img src="image_13.png" /><figcaption>**Figure 10.** 2D images of pedestal pressure gradient profiles of (a) H-mode and (b) super H-mode during the eELM crash. The transport caused by ELM in saturation phase leads to much more collapse at the pedestal region.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Contour plot of perturbed pressure from the beginning of saturation phase to the end of the pedestal collapse at the outer midplane in the H-mode. The horizontal axis is the normalized poloidal flux, and the vertical axis is the toroidal angle (1/5 torus). The collapse due to one ELM in H-mode is small, and the toroidal mode number of dominant mode transforms from low $n = 5$ to high $n = 10$.</figcaption></figure>

saturation time of ELM, and then it will determine the ELM size. It should be noted that $n = 0$ mode of $A_\|$ and $\phi$ are not considered in these simulations, which play important roles for the saturation of ELMs. Since the peeling mode is the current driven mode, without the $n = 0$ mode of $A_\|$, the energy loss level will be overestimated. In addition, if the mean flow is applied initially, then it contributes to the stabilization of the ballooning mode linearly. However, $n = 0$ mode of $\phi$ is generated to cancel out the mean flow in the nonlinear phase if it is included, so that the saturation level might be higher than that without the mean flow [45].

To describe the difference between the cases with small and large $E \times B$ velocity shear, a comparison of contour plots of perturbed pressure at the outer midplane during the saturation phase is illustrated in figure 18. As shown in figure 18, $n = 5$ mode is dominant in the case with $D_s = 10$, and $n = 10$ mode dominates in the case with $D_s = -10$ due to the less stabilizing effect on the high $n$ mode of smaller $E \times B$ velocity shear. The

instability saturates earlier with $D_s = -10$, and it causes less pedestal collapse than in the case with $D_s = 10$. It is worth noting that the physical process in the saturation phase plays an important role for us to understand the difference between the cases with $D_s = -10$ and $D_s = 10$. Based on the analysis of the nonlinear simulations in figure 18, it can be concluded that the smaller $E \times B$ velocity shear may be useful for mitigating the ELMs at the peeling boundary in super H-modes. However, the phenomenon with changes of $E \times B$ velocity shear may be different when the equilibria are located at the ballooning boundary.

### 3.3.2. Effect of $E \times B$ velocity shear on the ELM at ballooning boundary.
Larger $E \times B$ velocity shear may stabilize the ballooning mode and suppress the ELM. As shown by the green star in figure 19, the equilibrium with $n_{ped} \sim 8 \times 10^{19} m^{-3}$ that is located at the lower boundary is selected. The temperature
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Contour plot of perturbed pressure from the beginning of saturation phase to the end of the pedestal collapse at the outer midplane in super H-mode. The horizontal axis is the normalized poloidal flux, and the vertical axis is the toroidal angle (one-fifth of the torus). The collapse due to one ELM in super H-mode is very large, and the toroidal mode number of dominant mode does not change during the saturation phase.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 13.** (a) Linear growth rates of different $n$ and (b) evolutions of RMSP in the initial equilibrium ($D_s = 0$), equilibria with small shear ($D_s = -10$) and large shear ($D_s = 10$).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 14.** Evolutions of pedestal pressure profiles with (a) weak and (b) strong $E \times B$ velocity shear.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure 15.** Evolutions of (a) ELM size and (b) growth rate with different $E \times B$ velocity shear. Positive $D_s$ indicates the larger $E \times B$ velocity shear, and negative $D_s$ represents smaller $E \times B$ velocity shear.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 16.** Frequency of instabilities in saturation phase of the cases with (a) high $E \times B$ velocity shear ($D_s = 10$) and (b) low $E \times B$ velocity shear ($D_s = -10$). The frequency spectrum of the instabilities with high $E \times B$ velocity shear has a low value and only locates at one certificate region. But the frequency spectrum of the instabilities with low $E \times B$ velocity shear bifurcates into several parts.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 17.** The amplitudes of different toroidal number $n$ in the cases with (a) high $E \times B$ velocity shear ($D_s = 10$) and (b) low $E \times B$ velocity shear ($D_s = -10$). The amplitudes of each mode are normalized by standard deviation.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Figure 18.** Comparison of the contour plot of perturbed pressure at the outer midplane during the linear and saturation phase with small and large $E \times B$ velocity shear. $n = 10$ mode dominates in the case with $D_s = -10$, and $n = 5$ mode is dominant in the case with $D_s = 10$. The ELM size in the situation of large $E \times B$ velocity shear ($D_s = 10$) lead to much collapse at pedestal region.</figcaption></figure>

and bootstrap current in this case are lower than those in the equilibrium of super H-mode, which is illustrated above. In figure 20(a), it is obvious that the higher $n$ modes ($n \ge 15$) dominate and the larger $E \times B$ velocity shear can significantly reduce the growth rate of the higher $n$ modes. Then, with a very small growth rate (<0.03), the ELM, which is driven by ballooning modes, will have small amplitude ($\Delta W/W_{ped} < 3\%$) or just not occur (figure 20(b)). Moreover, if the $n = 0$ mode of $\phi$ is considered in this simulation, the RMS pressure may have a different tendency with the case located at the peeling boundary mentioned in section 3.3.1 [45].

For a clear understanding of the nonlinear simulation results at the ballooning boundary, the contour plot of the perturbed pressure is illustrated in figure 21. Figure 21 shows the comparison of contour plots of perturbed pressure in the outer midplane during the saturation phase with small and large $E \times B$ velocity shear. Based on this figure, we show that the ELM size is much larger with a small $E \times B$ velocity shear. The order of magnitude of the perturbed pressure is nearly 10–100 times higher in the case with small $E \times B$ velocity shear, and the collapse due to the ELM activity extends further inside the core of the plasma. Moreover, the pedestal profile

does not collapse with large $E \times B$ velocity shear due to the fact that the ballooning modes are well stabilized.

As mentioned in [4], the parameter space of the super H-mode is very narrow, and a single ELM crash may lead to a transition from super H-mode to H-mode, which is dominated by ballooning mode. Enlarging the parameter space of the super H-mode will be beneficial for sustaining the super H-mode operation. This study found that larger $E \times B$ velocity shear will stabilize the ELMs near the ballooning boundary and enlarge the parameter space of the super H-mode (figure 19). However, the increase in the ELM size due to the large $E \times B$ velocity shear at the peeling boundary may cause more energy loss due to a single ELM crash, which will make the transition from super H-mode to H-mode (ballooning mode dominates) easier. A summary is made of the nonlinear simulations with various equilibria and $E \times B$ velocity shear, which are shown in table 2.

By non-inductive means, for example, NBI [46, 47] and electron cyclotron resonance heating [48], the $E \times B$ velocity shear can be modulated in both negative and positive directions. Moreover, impurity seeding techniques such as laser blow off [49, 50] and super molecular beam injection [51] also
<figure><img src="image_22.png" /><figcaption>**Figure 19.** Parameter space of the HL-3 equilibrium, which is dominated by ballooning mode. ELM size may not be mitigated by larger $E \times B$ velocity shear at the peeling boundary, whereas the larger $E \times B$ velocity shear can suppress the ELM at ballooning boundary.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 20.** (a) The amplitude of different toroidal number $n$ in the equilibria near the ballooning boundary with low $E \times B$ velocity shear ($D_s = -10$) and high $E \times B$ velocity shear ($D_s = 10$). (b) Evolutions of RMSP with different $E \times B$ velocity shear. $\gamma/\omega_A = 0$ means the growth rate is close to or less than 0. The RMSP with large shear is much smaller than others.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Figure 21.** Comparison of contour plot of perturbed pressure at the outer midplane during the saturation phase with small and large $E \times B$ velocity shear at ballooning boundary. It should be noted that the amplitude of perturbed pressure with large $E \times B$ velocity shear is 10–100 times smaller than that with small $E \times B$ velocity shear.</figcaption></figure>
<table><caption>**Table 2.** Summary of the nonlinear simulation results with low/high *E* × *B* velocity shear at peeling/ballooning boundary.</caption><thead><tr><th rowspan="2">Location</th><th colspan="3">Low E × B shear</th><th colspan="3">High E × B shear</th></tr><tr><th>Dominant mode</th><th>ELM size</th><th>Growth rate</th><th>Dominant mode</th><th>ELM size</th><th>Growth rate</th></tr></thead><tbody><tr><td>Peeling boundary</td><td>Low n = 5</td><td>Small</td><td>High</td><td>Low n = 5</td><td>Large</td><td>High</td></tr><tr><td>Ballooning boundary</td><td>High n = 10</td><td>Large</td><td>High</td><td>Suppression</td><td>Suppression</td><td>Suppression</td></tr></tbody></table>

have the ability to change the *E* × *B* velocity shear. According to table 2, the ELM control method can be modified in different parameter spaces.

## 4. Summary

Based on the scan results by EPED, the HL-3 super H-mode regime may operate under the condition of $n_{ped} > 7 \times 10^{19} \text{ m}^{-3}$ (0.78 $n_{GW}$). Larger $B_T$ leads to much higher pedestal pressure, and the threshold of $n_{ped}$ for achieving the super H-mode regime is also increased. An ELM crash in the HL-3 super H-mode operation will lead to large energy loss in a very short time. It is essential to develop the ELM control method for super H-mode operation.

Using nonlinear simulation results of $E_r$ shear and edge current density, we reveal that both the *E* × *B* velocity shear and edge current density affect the growth rate and the time of linear phase $t_l$. The combination of the growth rate and $t_l$ determines the ELM size. Larger *E* × *B* velocity shear leads to smaller growth rates but longer $t_l$ with the equilibrium of the super H-mode regime near the peeling boundary. A lower current density decreases both of these two parameters and mitigates the ELM size. The development of a technique that can decrease the current and maintain the pedestal pressure is needed. Moreover, enlarging the parameter space of super H-mode will be beneficial for sustaining the super H-mode operation. This study found that a larger *E* × *B* velocity shear will stabilize the ELMs near the ballooning boundary and enlarge the parameter space of the super H-mode. However, the increase in the ELM size due to the large *E* × *B* velocity shear at the peeling boundary may cause more energy loss due to a single ELM crash, which will make the transition from super H-mode to H-mode (ballooning mode dominates) easier.

According to the simulation results of this paper, it looks promising that HL-3 achieves the super H-mode in experiments. Furthermore, HL-3 will try to explore more physical details of the super H-mode regime and provide references to high-performance operations for future devices such as ITER.

## Acknowledgments

The author wishes to thank for the guidance by V.C. Chan. We would like to express our thanks to the permission to use the codes, EFIT, ELITE, EPED, ONETWO and OMFIT framework from the General Atomics theory group, BOUT ++

developed by LLNL and the University of York, and the guidance of EUROPED by S. Saarelma. This work is mainly supported by the National Natural Science Foundation of China under Grant Nos. 11922503 and U21A20440. It is also partially supported by the National Key R&D Program of China under Grant No. 2017YFE0301106.

## ORCID iDs

Yunpeng Zou https://orcid.org/0000-0003-3311-5931

Min Xu https://orcid.org/0000-0001-3059-7026

## References

[1] Aymar R., Barabaschi P. and Shimomura Y. 2002 The ITER design *Plasma Phys. Control. Fusion* **44** 519

[2] Wan Y. et al 2017 Overview of the present progress and activities on the CFETR *Nucl. Fusion* **57** 102009

[3] Doyle E.J. et al 2007 Chapter 2: Plasma confinement and transport *Nucl. Fusion* **47** S18

[4] Snyder P.B. et al 2015 Super H-mode: theoretical prediction and initial observations of a new high performance regime for tokamak operation *Nucl. Fusion* **55** 083026

[5] Snyder P.B. et al 2019 High fusion performance in super H-mode experiments on alcator C-mod and DIII-D *Nucl. Fusion* **59** 086017

[6] Knolker M. et al 2020 Optimizing the super H-mode pedestal to improve performance and facilitate divertor integration *Phys. Plasmas* **27** 102506

[7] Ding S., Garofalo A.M., Knolker M., Marinoni A., McClenaghan J. and Grierson B.A. 2020 On the very high energy confinement observed in super H-mode DIII-D experiments *Nucl. Fusion* **60** 034001

[8] Soloman W.M. et al 2016 Exploration of the Super H-mode regime on DIII-D and potential advantages for burning plasma devices *Phys. Plasmas* **23** 056105

[9] Ernst D.R. et al 2000 Transitionless enhanced confinement and the role of radial electric field shear *Phys. Plasmas* **7** 615

[10] Manini A., Angioni C., Peeters A.G., Ryter F., Jacchia A., Maggi C.F. and Suttrop W. (the ASDEX Upgrade Team) 2006 Role of $T_e/T_i$ and $\nabla v_{tor}$ in ion heat transport of ASDEX Upgrade H-mode plasmas *Nucl. Fusion* **46** 1047

[11] Lang P.T. et al 2013 ELM control strategies and tools: status and potential for ITER *Nucl. Fusion* **53** 043004

[12] Liang Y. et al 2015 Edge localized mode *Active Control of Magnetohydrodynamic Instabilities in Hot Plasmas* (Springer Series on Atomic, Optical and Plasma Physics) vol 83 (Springer) ch 5 (available at: www.springer.com/series/0411)

[13] Viezzer E. et al 2018 Access and sustainment of naturally ELM-free and small-ELM regimes *Nucl. Fusion* **58** 115002
[14] Zeyu L. et al 2021 Edge localized mode characteristics and divertor heat flux during stationary and transient phase for CFETR hybrid scenario *Plasma Phys. Control. Fusion* **63** 035006

[15] Xue L. et al 2020 Integrated plasma scenario analysis for the HL-2M tokamak *Nucl. Fusion* **60** 016016

[16] Duan X.R. et al 2022 Progress of HL-2A experiments and HL-2M program *Nucl. Fusion* **62** 042020

[17] Snyder P.B., Groebner R.J., Leonard A.W., Osborne T.H. and Wilson H.R. 2009 Development and validation of a predictive model for the pedestal height *Phys. Plasmas* **16** 056118

[18] Saarelma S., Challis C.D., Garzotti L., Frassinetti L., Maggi C.F., Romanelli M. and Stokes C. 2018 Integrated modelling of H-mode pedestal and confinement in JET-ILW *Plasma Phys. Control. Fusion* **60** 014042

[19] Li K. et al 2020 Study of H-mode pedestal predictive model on EAST tokamak *Plasma Phys. Control. Fusion* **62** 115007

[20] Li K., Lao L.L., Li G.Q., McClenaghan J., Jian X., Osborne T., Evans T.E., Smith S.P., Meneghini O.M. and Snyder P.B. 2021 Study of H-mode pedestal model for helium plasmas in DIII-D *Nucl. Fusion* **61** 096002

[21] Snyder P.B., Groebner R.J., Hughes J.W., Osborne T.H., Beurskens M., Leonard A.W., Wilson H.R. and Xu X.Q. 2011 A first-principles predictive model of the pedestal height and width: development, testing and ITER optimization with the EPED model *Nucl. Fusion* **51** 103016

[22] Snyder P.B., Wilson H.R., Ferron J.R., Lao L.L., Leonard A.W., Mossessian D., Murakami M., Osborne T.H., Turnbull A.D. and Xu X.Q. 2004 ELMs and constraints on the H-mode pedestal: peeling-ballooning stability calculation and comparison with experiment *Nucl. Fusion* **44** 320–8

[23] Doyle E.J. et al 2010 Demonstration of ITER operational scenarios on DIII-D *Nucl. Fusion* **50** 075005

[24] Wilson H.R. and Miller R.L. 1999 Access to second stability region for coupled peeling-ballooning modes in tokamaks *Phys. Plasmas* **6** 3

[25] Zhu Y.-R., Li Z.-Y., Chan V.S., Chen J.-L., Jian X., Dudson B.D., Garofalo A.M., Snyder P.B., Xu X.-Q. and Zhuang G. 2020 Achieving a robust grassy-ELM operation regime in CFETR *Nucl. Fusion* **60** 046014

[26] Meneghini O. et al 2016 Integrated fusion simulation with self-consistent core-pedestal coupling *Phys. Plasmas* **23** 042507

[27] Wang X., Mordijck S., Doyle E.J., Zeng L., Staebler G.M., Meneghini O. and Smith S.P. 2018 Role of turbulence regime on determining the local density gradient *Nucl. Fusion* **58** 016025

[28] Kang J. et al 2020 Role of fast-ion transport manipulating safety factor profile in KSTAR early diverting discharges *Nucl. Fusion* **60** 126023

[29] Kinsey J.E., Lao L.L., Meneghini O., Candy J., Snyder P.B. and Staebler G.M. 2022 Equilibrium reconstruction of DIII-D plasmas using predictive modeling of the pressure profile *Phys. Plasmas* **29** 062502

[30] Gong X. et al 2019 Integrated operation of steady-state long-pulse H-mode in experimental advanced superconducting tokamak *Nucl. Fusion* **59** 086030

[31] Xu X.Q., Dudson B.D., Snyder P.B., Umansky M.V., Wilson H.R. and Casper T. 2011 Nonlinear ELM simulations based on a nonideal peeling-ballooning model using the BOUT++ code *Nucl. Fusion* **51** 103040

[32] Xu X.Q., Dudson B., Snyder P.B., Umansky M.V. and Wilson H. 2010 Nonlinear simulations of peeling-ballooning modes with anomalous electron viscosity and their role in edge localized mode crashes *Phys. Rev. Lett.* **105** 175005

[33] Seto H., Xu X.Q., Dudson B.D. and Yagi M. 2024 Two-stage crash process in resistive drift ballooning mode driven ELM crash *Phys. Plasmas* **31** 03251

[34] Ebrahimi F. and Bhattacharjee A. 2023 Plasmoid-mediated reconnection during nonlinear peeling-ballooning edge-localized modes *Nucl. Fusion* **63** 126042

[35] Pamela S., Bhole A., Huijsmans G.T.A., Nkonga B., Hoelzl M., Krebs I. and Strumberger E. 2020 Extended full-MHD simulation of non-linear instabilities in tokamak plasmas *Phys. Plasmas* **27** 102510

[36] Li Z.-Y. et al 2018 Ideal MHD stability and characteristics of edge localized modes on CFETR *Nucl. Fusion* **58** 016018

[37] Polevoi A.R., Loarte A., Dux R., Eich T., Fable E., Coster D., Maruyama S., Medvedev S.Y., Köchl F. and Zhogolev V.E. 2018 Integrated simulations of H-mode operation in ITER including core fuelling, divertor detachment and ELM control *Nucl. Fusion* **58** 056020

[38] Waelbroeck F.L. and Chen L. 1991 Ballooning instabilities in tokamaks with sheared toroidal flows *Phys. Fluids B* **3** 601

[39] Miller R.L., Waelbroeck F.L., Hassam A.B. and Waltz R.E. 1995 Stabilization of ballooning modes with sheared toroidal rotation *Phys. Plasmas* **2** 3676

[40] Terry P.W. 2000 Suppression of turbulence and transport by sheared flow *Rev. Mod. Phys.* **72** 109–65

[41] Simakov A.N. and Catto P.J. 2003 Drift-ordered fluid equations for field-aligned modes in low-β collisional plasma with equilibrium pressure pedestals *Phys. Plasmas* **10** 4744–57

[42] Xi P.W., Xu X.Q., Wang X.G. and Xia T.Y. 2012 Influence of equilibrium shear flow on peeling-ballooning instability and edge localized mode crash *Phys. Plasmas* **19** 092503

[43] Zhu Y. et al 2022 Simulations on edge localized modes mitigation with impurity seeding in the HL-2A tokamak *Nucl. Fusion* **62** 076011

[44] Seto H., Dudson B.D., Xu X.-Q. and Yagi M. 2023 A BOUT++ extension for full annular tokamak edge MHD and turbulence simulations *Comput. Phys. Commun.* **283** 108568

[45] Seto H., Xu X.Q., Dudson B.D. and Yagi M. 2020 Impact of equilibrium radial electric field on energy loss process after pedestal collapse *Contrib. Plasma Phys.* **60** 201900158

[46] Maggi C.F. et al 2011 L-H threshold at low density and low momentum input in the JET tokamak 38th EPS Conf. on *Plasma Physics* (Strasbourg Convention Centre, France, 27 June–1 July 2011) vol P4 p 061

[47] Chen X. et al 2020 Expanding the parameter space of the wide-pedestal QH-mode towards ITER conditions *Nucl. Fusion* **60** 092006

[48] Reisner M. et al 2021 Investigation of increased core ion temperatures in high-beta advanced scenarios in ASDEX Upgrade 47th EPS Conf. on *Plasma Physics* (Online, 21–25 June 2021) vol O5 p 102 (available at: https://epsplasma2020.eu/)

[49] Zhang Y.P. et al 2018 Control of edge localized modes by pedestal deposited impurity in the HL-2A tokamak *Nucl. Fusion* **58** 046018

[50] Xiao G.L. et al 2021 Dual effects of impurity seeding on pedestal turbulence and ELMs in the HL-2A tokamak *Nucl. Fusion* **61** 116011

[51] Zhong W.L. et al 2019 Impact of impurity mixture gas seeded by supersonic molecular beam injection on edge-localized modes in the HL-2A tokamak *Nucl. Fusion* **59** 076033