Resistive reduced MHD modeling of multi-edge-localized-mode

The full dynamics of a multi-edge-localized-mode (ELM) cycle is modeled for the first time in realistic tokamak X -point geometry with the nonlinear reduced MHD code JOREK . The diamagnetic rotation is found to be instrumental to stabilize the plasma after an ELM crash and to model the cyclic reconstruction and collapse of the plasma pressure profile. ELM relaxations are cyclically initiated each time the pedestal gradient crosses a triggering threshold. Diamagnetic drifts are also found to yield a near-symmetric ELM power deposition on the inner and outer divertor target plates, consistent with experimental measurements

Introduction.-Edge localized modes (ELMs) are magnetohydrodynamic (MHD) instabilities occurring in fusion plasmas in high confinement regime (H-mode) of tokamaks and occasionally in stellerators, and are sometimes compared in the early stage of their development with the solar flares occurring in astrophysical plasmas [1]. As the ELMs represent a particular concern for the lifetime of the divertor in ITER [2], their dynamics has been widely studied experimentally and theoretically [3][4][5][6]. In X-point tokamak plasmas, MHD instabilities develop due to a large edge pressure gradient (so-called ballooning modes) or a large edge current (so-called peeling-tearing modes), triggering the formation of ELM filaments. The ELMs consist in the quasiperiodic ejection of these filaments in H-mode regimes across the transport barrier at the plasma edge (also called pedestal).
So far, nonlinear MHD computations with codes such as M3D [7], BOUT++ [8,9], NIMROD [10], and JOREK [11][12][13] have been able to reproduce a single ELM crash in simulations, in realistic tokamak geometry. Other works using a model including two-fluid diamagnetic effects [14,15] have shown transport bursts comparable to ELM relaxations in simplified cylindrical geometry. This Letter presents the first simulations of a multi-ELM cycle in realistic toroidal geometry, using the nonlinear reduced MHD code JOREK [11] recently extended to add the two-fluid diamagnetic effects [16]. The diamagnetic effects are known to have a stabilizing effect on plasma instabilities. They are shown to prevent the magnetic field from reconnecting and/or to reduce the growth rate of ideal and resistive instabilities in both astrophysical [17] and laboratory plasmas [18][19][20]. In particular, in tokamaks, they are evidenced to be a key parameter for simulating cycles of sawtooth crashes [21]. We show here that the diamagnetic rotation is also a key parameter for modeling the ELM cycling dynamics: as the diamagnetic effects stabilize peeling-ballooning instabilities after an ELM crash, they allow for the pedestal to rebuild, thus leading to an ELM cycle. This approach is different and complementary to that in Ref. [22] in which it is an externally imposed E × B stabilization in the plasma edge that allows for oscillation relaxations.
Multi-ELM cycle simulations.-Experimentally, ELM relaxations consist in a cyclical phenomenon, resulting in the quasiperiodic deposition of energy on the divertor target plates. Modeling this phenomenon requires to reproduce this cyclical behavior. Besides, the simulation of ELM cycles rather than a single ELM crash, which depends on the initial unstable pressure profile, involves different physical mechanisms. Indeed, after the first ELM crash, memory of the choice for the initial pressure profile is lost. The phase coherence between modes that determines the ELM instability growth dynamically evolves, nonlinearly, through the heat, particle, and current sources. And so is the energy distribution among the toroidal harmonics. In this respect, an ELM crash starting from such a coherent state is significantly different from a first initial relaxation triggered by the choice of an initial state. This Letter discusses the specificities of the multiple ELM dynamics starting from consistent inter-ELM states with respect to the usual single ELM relaxation from an initial chosen condition.
As we consider long runs of ELM-cycle simulations (>10 4 Alfvén times ∼ 10 −2 s), we have to deal with very different time scales: the diffusive time scale τ D ∼ a 2 =D ⊥ ∼10 −1 − 1 s (D ⊥ being the perpendicular diffusion coefficient), the time scale of the pedestal reconstruction after an ELM (∼ 10 −3 − 10 −1 s), and the ELM crash time scale (∼10 −4 s). Thus, the presence of heat and particle sources-needed to rebuild the pedestal after an ELM-as much as diffusion propagating the heat and particles, is essential in the simulations. The heat source is chosen to follow the initial temperature profile and the particle source is taken constant over the plasma. Heat and particle diffusion coefficients are reduced in the pedestal in order to preserve a transport barrier at the plasma edge and a 3-cm-large pedestal. The injected heat and particle sources determine the pedestal reconstruction time scale and thus, the simulated ELM frequency. As for the current profile, it is driven back close to the initial realistic current profile J 0 (including Ohmic and bootstrap components) after an ELM crash, using the Krook operator ηðJ − J 0 Þ in the Ohms law (η being the plasma resistivity). An ongoing work aims at implementing a self-consistently evolving bootstrap current.
In this Letter, ELM cycles are modeled for Joint European Torus (JET)-like plasma parameters and geometry, similar to Ref. [16]: major radius R 0 ¼ 3 m, minor radius a ¼ 1 m, toroidal magnetic field B t ¼ 2.9 T, and safety factor q 95 ∼ 3. Established H-mode experimental profiles are taken initially, with central electron density n e;0 ¼ 6 × 10 19 m −3 and central temperature T e;0 ¼ 5 keV. The pedestal density and temperature are n e;ped ¼ 3.8 × 10 19 m −3 and T e;ped ¼ 2.5 keV. Ion and electron temperatures are assumed to be equal. Because of computational restrictions, the central resistivity is taken η 0 ¼ 2.5 × 10 −7 Ω · m, which is 2 orders of magnitude larger than the Spitzer resistivity η Spitzer ¼ 2.5 × 10 −9 Ω · m. The resistivity profile follows a T −3=2 radial dependence. A JETrealistic diamagnetic velocity is taken for both species s (ions and electrons):Ṽ Ã s ¼ Z s ðR=R 0 Þ 2 τ IC =ρ·b × ∇P s . Z s and P s are, respectively, the charge number and the scalar pressure of the species s, R is the horizontal coordinate, ρ is the mass density, andb is the magnetic field normalized to B t . The realistic value of the diamagnetic parameter (inverse of the ion cyclotron frequency) τ IC ¼ m i =eB t ∼ 7 × 10 −9 s is taken. The equilibrium value of the diamagnetic velocities in the pedestal is m=s, which is in the ballpark of the experimental observations [23]. The full description of the model is given on Ref. [16]. Compared to Ref. [16], the source of parallel rotation S V jj and the neoclassical friction (∇ ·Π i;neo term) are not included in the modeling to focus only on the impact of the diamagnetic flow on the ELM dynamics.
A simulation is run as follows: once the axisymmetric (n ¼ 0) equilibrium flows are established (mainly due to diamagnetic drifts and sheath conditions [16]), the toroidal modes of perturbation n ¼ 2; 4; 6, and 8 are added in simulation. Linear runs show that the n ¼ 8 mode is the most unstable mode in this configuration, as the larger n > 9 modes are stabilized by the diamagnetic effects. The grid resolution (102 points in the radial direction, 128 in the poloidal one) is chosen such that the mode structures up to n ¼ 8 are well resolved. Further coupling with high n > 9 modes is likely to happen as well as coupling with modes of odd parity and should be added if more quantitative predictions are required. Quantitatively addressing these questions during an ELM cycle is, however, currently beyond the scope of our numerical capabilities. In the present Letter, the focus is on the ELM behavior in quasiperiodic cycles. Interestingly, it is found to significantly differ from the oft-computed single-ELM crash.
In ELM simulations without diamagnetic effects, the ballooning and/or peeling-tearing instabilities at the onset on the ELM crash generate the stochastization of the edge magnetic field. Yet after the ELM crash, the MHD activity is not fully stabilized and the edge magnetic field remains stochastic: the residual MHD activity still generates an enhanced edge transport that does not allow for the pedestal reconstruction. Therefore, a second ELM crash cannot be obtained as the plasma remains below the peeling-ballooning stability limit. An example is given for a single n ¼ 8 ELM (only the n ¼ 0 and n ¼ 8 modes are included in the simulation) modeled without taking into account the diamagnetic effects ( Fig. 1): during the ELM crash, characterized by a large peak of kinetic and magnetic energy of the n ¼ 8 mode (Fig. 1), the plasma edge is stochastized in a large extent: the region for normalized flux ψ norm > 0.85 is fully reconnected [Figs. 2(a) and 2(c)]. After the crash, the instability does not completely vanish (Fig. 1) and the magnetic field remains partially reconnected: large magnetic islands subsist for ψ norm > 0.85 and an ergodic layer remains for ψ norm > 0.95 [Figs. 2(b) and 2(d), given for t∼ 7 × 10 3 t A ]. Thus, the edge transport remains high due to both the stochastic layer enhancing the heat parallel diffusivity and the E × B convection of particles. However, the addition of the diamagnetic rotation in modeling has two major effects: first, the amplitude of the ELM crash is smaller due to the diamagnetic stabilization, and second, the n ¼ 8 mode is completely damped after the crash. The plasma is actually stabilized by the diamagnetic effects after the crash, but is not destabilized again until the pedestal pressure gradient that increases due to the injected heat source reaches the ELM-triggering threshold. Then, a second ELM crash occurs, and so forth. The simulation of a multi-ELM cycle (multiharmonic simulation n ¼ 0 to 8 by step 2) is presented in Fig. 3. The ELM energy is one order of magnitude smaller than in Fig. 1 due to the diamagnetic stabilization. The ELM cycle can be decomposed into two parts: transient ELMs-from t ¼ 3200t A to t ¼ 6200t Aare followed by a quasiperiodic repetitive regime. The system evolves from a state depending on the initial conditions towards a quasiperiodic cycling regime which mainly depends on the injected heat and particles and on the diamagnetic frequency. The three first transient ELMs are largely dominated by the most unstable mode (first, n ¼ 8 and then, n ¼ 6) and there is only a few coupling between modes. However, after four ELMs, we get a quasiperiodic regime driven by the n ¼ 2-8 modes strongly nonlinearly coupled with each other. Each ELM in this regime is initiated by the growth of the n ¼ 6 mode directly followed by the coupled growth of the other modes.
After a crash, because of the quick refilling of the pedestal by the heat and particle sources, the plasma is barely stabilized by the diamagnetic effects before being destabilized again by the increase in pressure gradient: this leads to a large frequency (f ELM ∼ 3 kHz, f ELM normalized to Afvén time ∼ 1.6 × 10 −3 ) of small ELMs. Each ELM in the quasiperiodic regime is characterized by the collapse of the edge pressure gradient, presented in Fig. 4. This pedestal relaxation occurs when a threshold in pressure gradientcorresponding to the ELM triggering threshold-is reached. The peeling-ballooning diagram (maximal pressure gradient as a function of the mean pedestal current Fig. 5) shows that these ELM are constituted of three steps: the linear growth of the instability as the pressure gradient increases, then the ELM crash-occurring at approximately the same threshold for all the ELMs in quasiperiodic regime-and last, the relaxation of the density and temperature profiles. The increase of the pressure gradient  in the pedestal due to the heat source is going together with the increase of the rotational shear at the plasma edge, due to the increase in diamagnetic velocity proportional to the pressure gradient. At the ELM crash, an n ¼ 0 flow is induced by Maxwell stress [24], increasing the edge rotational shear. The plasma filaments are then sheared by this (n ¼ 0) flow and expelled out of the plasma. Near-symmetric deposition on divertor.
-After a small delay (∼ 100t A ) corresponding to the ion parallel time τ jj;i ¼ πqR=c s (where c s is the sound speed), the energy and particles due to the ELM filaments reach the divertor target plates at the sound speed. The particles are convected by E × B and diamagnetic drifts while the energy is mainly transported by parallel conduction. Although the power deposition due to the four first ELMs is very variable, the integrated peak power deposited by an ELM on the divertor (Fig. 6) is approximately the same for all the ELMs in quasiperiodic regime (∼ 5-6 MW on the outer divertor, ∼2-3 MW on the inner divertor). This again shows that the first ELMs before reorganization are singular. Note that in quasiperiodic regime, the peak power deposition on the inner divertor is slightly delayed (by ∼ 50t A ) compared to the outer deposition: because of the ballooning character of the instabilities, the filaments are mostly expelled from the low field side, so the filaments first hit the outer target.
Furthermore, the deposition on the inner and outer divertor targets is observed to be near symmetric in this simulation: the integrated deposited power has the same order of magnitude on both targets, even though roughly 2 times more power is deposited on the outer divertor target. Experimentally, the peak power deposition is either symmetric on inner or outer targets, or the inner target receives twice more power [25,26] after a crash. In previous modeling performed without diamagnetic drifts, the outer target received almost all the ELM power deposition, which was contrary to the experimental observations. An n ¼ 2-8 ELM simulation is performed with realistic JET parameters corresponding to the shot #77329 described in [27] for two different cases: one without including flows in the model (full line in Fig. 7) and one with diamagnetic effects, neoclassical friction and toroidal source of rotation included (dash line in Fig. 7). In the case without flows, almost all the heat flux generated by the ELM filaments is deposited on the outer divertor, whereas in the case where flows are included in the model, the deposited heat flux is near symmetric on the inner and outer divertor targets (Fig. 7). Therefore, even though the filaments are mostly expelled from the outer region and the heat transport is mainly diffusive in direction of the outer region, the particle convection towards the inner region is enhanced by the diamagnetic flow, leading to a near-equal repartition of the power (proportional to density and temperature) on inner and outer divertor plates. Hence, plasma flows and, in particular, the two-fluid diamagnetic drifts allow for simulating more realistically the ELM cyclical dynamics up to the deposition on divertor.   Discussion and conclusions.-Even though computational restrictions impose to choose in the simulation a plasma resistivity that is 2 orders of magnitude larger than the realistic value, the modeling of these resistive ELMs represents a significant step forward in the understanding of the cycling dynamics. The addition of the two-fluid diamagnetic rotation in MHD models allows for the stabilization of the plasma after an ELM relaxation. The plasma reorganizes and memory of the initial pressure profile is lost. The phasing and the energy repartition between modes is consistently determined. The steepening of the pressure profile generated by the pedestal reconstruction destabilizes again the edge plasma until the edge pressure gradient reaches the ELM-triggering threshold: a new ELM crash then occurs. Similar coupling between modes, similar maximal pressure gradient reached when the crash occurs, and similar power deposition on the divertor plates are cyclically recovered for all the ELMs in the quasiperiodic regime. These differ much from the first transient ELMs, pointing out the importance of simulating cycles rather than a single ELM crash. In addition, the inclusion of the diamagnetic drifts induces a nearsymmetric ELM power deposition on the inner and outer divertor target plates, which is in closer agreement with experimental measurements. This work has been carried out within the framework of the EUROfusion Consortium. It has received funding from the European Union Horizon 2020 research and innovation programme under grant agreement number 633053 and from the National French Research Program (ANR): ANEMOS (2011) and E2T2 (2010). This work was granted access to the HPC resources of Aix-Marseille Université financed by the project Equip@Meso (ANR-10-EQPX-29-01). A part of this work was carried out using the HELIOS supercomputer system (IFERC-CSC), Aomori, Japan, under the Broader Approach collaboration, implemented by Fusion for Energy and JAEA, and using the CURIE supercomputer, operated into the TGCC by CEA, France, in the framework of GENCI and PRACE projects. The views and opinions expressed herein do not necessarily reflect those of the European Commission or the ITER Organization.