Recent progress in the quantitative validation of JOREK simulations of ELMs in JET

This content has been downloaded from IOPscience. Please scroll down to see the full text. Download details: IP Address: 207.162.240.147 This content was downloaded on 16/05/2017 at 13:36 Please note that terms and conditions apply. Recent progress in the quantitative validation of JOREK simulations of ELMs in JET View the table of contents for this issue, or go to the journal homepage for more 2017 Nucl. Fusion 57 076006 (http://iopscience.iop.org/0029-5515/57/7/076006) Home Search Collections Journals About Contact us My IOPscience 1 © 2017 Culham Centre for Fusion Energy Printed in the UK 1. Introduction and motivation Elaborate experimental scaling laws for the ITER tungsten divertor provide an estimate of the ELM (Edge–Localised– Mode) and inter-ELM target heat fluxes [1, 2]. These predictions can be reinforced by numerical simulations of large-scale instabilities, like peeling–ballooning (PB) modes, to describe the characteristic dynamics of ELMs, as well as small-scale turbulence of kinetic-ballooning modes (KBMs) and ion-temperature-gradient (ITG) instabilities to describe the cross-field transport in the pedestal which regulates the plasma and energy exhaust across the separatrix in inter-ELM regimes. Several nonlinear MHD codes, such as JOREK, BOUT++, HESEL and EMEDGE3D in Europe [3–8], M3D-C1 and NIMROD in the US [9, 10], or MEGA in Japan [11–13], can obtain advanced ELM simulations, with challenging physics effects like bi-fluid diamagnetic rotation and current, low resistivity and viscosity levels, as well as high poloidal/ toroidal resolutions. In the last decade, nonlinear MHD codes Nuclear Fusion Recent progress in the quantitative validation of JOREK simulations of ELMs inJET S.J.P.Pamela1, G.T.A.Huijsmans2,3, T.Eich4, S.Saarelma1, I.Lupelli1, C.F.Maggi1, C.Giroud1, I.T.Chapman1, S.F.Smith1,5, L.Frassinetti6, M.Becoulet2, M.Hoelzl4, F.Orain4, S.Futatani7 and JET Contributorsa EUROfusion Consortium, JET, Culham Science Centre, Abingdon, OX14 3DB, United Kingdom 1 CCFE, Culham Science Centre, Abingdon, Oxon, OX14 3DB, United Kingdom 2 CEA, IRFM, F-13108 Saint-Paul-lez-Durance, France 3 Eindhoven University of Technology, Eindhoven, Netherlands 4 Max-Planck-Institut für Plasmaphysik, Boltzmannstrasse 2, D-85748 Garching bei München, Germany 5 Department of Physics, York Plasma Institute, University of York, Heslington, York YO10 5DD, United Kingdom 6 Division of Fusion Plasma Physics, KTH Royal Institute of Technology, Stockholm SE, Sweden 7 CASE Department, Barcelona Supercomputing Center, Barcelona, Spain E-mail: [email protected] Received 19 December 2016, revised 31 March 2017 Accepted for publication 20 April 2017 Published 16 May 2017 Abstract Future devices like JT-60SA, ITER and DEMO require quantitative predictions of pedestal density and temperature levels, as well as inter-ELM and ELM divertor heat fluxes, in order to improve global confinement capabilities while preventing divertor erosion/melting in the planning of future experiments. Such predictions can be obtained from dedicated pedestal models like EPED, and from non-linear MHD codes like JOREK, for which systematic validation against current experiments is necessary. In this paper, we show progress in the quantitative validation of the JOREK code using JET simulations. Results analyse the impact of diamagnetic terms on the dynamics and size of the ELMs, and evidence is provided that the onset of type-I ELMs is not governed by linear MHD stability alone, but that a nonlinear threshold could be responsible for large MHD events at the plasma edge. Keywords: JET-ILW, ELM, pedestal, MHD, peeling ballooning modes, JOREK, simulations (Some figuresmay appear in colour only in the online journal) S.J.P. Pamela etal Recent progress in the quantitative validation of JOREK simulations of ELMs in JET Printed in the UK 076006 NUFUAU © 2017 Culham Centre for Fusion Energy 57 Nucl. Fusion NF 10.1088/1741-4326/aa6e2a Paper 7 Nuclear Fusion IOP International Atomic Energy Agency a See author list of Litaudon [50] 2017 1741-4326 1741-4326/17/076006+11$33.00 https://doi.org/10.1088/1741-4326/aa6e2a Nucl. Fusion 57 (2017) 076006 (11pp) S.J.P. Pamela etal 2 have focused their efforts on obtaining qualitative agreement with experimental observations, by considering key ELM characteristics like the formation of hot plasma filaments that are ejected through the separatrix, the collapse of the pedestal pressure, and the transport of energy to the divertor and first wall components. Recent increases in computational resources available to fusion research are now enabling the start of quanti tative comparisons with experiments, which is the compulsory path towards predictions for future devices like ITER. In sight of the urgent need for such predictions, this paper discusses the latest progress of quantitative validations for the nonlinear MHD code JOREK against experimental data from the JET-ILW device. In particular, one of the main issues for all nonlinear MHD codes that address ELM physics is the methodology used for the simulation set-up, and more precisely the initial boundary conditions. Unless multiple ELM cycles are simulated, the initial condition is typically chosen to be an unstable axisymmetric equilibrium with small toroidal perturbations (of the order of numerical noise). This paper will show that, if these initial conditions are used, nonlinear coupling cannot be obtained in the critical phase of the ELM crash; it will only occur in the late phase of the ELM, at which point it is too late to make a significant difference in ELM dynamics. The paper is organised as follows. In section2, we provide details on the JOREK code, including the physics model used for the simulations, and how JET experimental data is used to create initial conditions for the simulations. In section3, we describe the results on the quantitative validation of JOREK simulations for ELMs in JET, and how these are affected by diamagnetic terms. Section3 will also discuss the nonlinear aspect of simulations and the role played by initial value conditions concerning coupling between toroidal modes. In the conclusion, section4, we discuss how predictions for future devices with the JOREK code need to rely on multi-cylce ELM simulations, and what this implies for future simulations. 2. The JOREK code 2.1. The physics and numerical models The 3D nonlinear MHD code JOREK was developed by Huysmans et al with the specific aim to produce simulations of edge-localised-modes [3, 4]. The MHD model used for the present paper is similar to that used in previous ELM studies [14–18]. It is a five-field reduced MHD model for the variables ψ (poloidal magnetic flux), Φ (electric potential), v∥ → (parallel velocity), ρ (density), T (total temperature), including the two-fluid diamagnetic effects. The reduction of the equationsassumes that the perpendicular velocity lies in the poloidal plane, and that the toroidal magnetic field is constant in time, so that the total plasma velocity and the total magnetic field are expressed respectively as vvvvvv ∥∥ ∥ →→→→→→ → →→ δ ρ =+=++ =+×∇Φ+ ×∇ φφ ⊥∗ ∗ vB Re R ep , Ei i tot (1) →→ → →→ ψ=+=+∇× φφ φ BBB F R e R e 1, po (2) where R is the major radius, → φ e is the toroidal unit vector and =FBR ooo , with Bo being the magnetic field amplitude at the reference major radius R = Ro. The diamagnetic component of the perpendicular velocity is represented by the third term v →→ δρ=×∇ φ∗ ∗− Re p i i 1 , where pi is the ion pressure and ()δ=Ω ∗− R ci o1 , with the ion gyrofrequency /Ω=eB m ci oi . Substituting the identities (1) and (2) into the visco-resistive MHD equationsgives the reduced MHD model, first derived by Strauss [19], with two equationsfor the parallel and the perpendicular momentum, v vv vv vv () () → →→ →→ →→ →→ ρρ µµ =− ⋅∇ −∇ +× +∇ ++ ∇+ ∗⊥ ∗∗ tpJB d d , EiE Ei Ei 2hyp 4 (3) v vv vv ∥ ∥∥∥∥ ∥ → →→ →→ ρρ µµ=− ⋅∇ −∇ +∇ +∇ t p d d, 2hyp 4 (4) ()[] [] ψ ηψ φ δ ρψδ ρφ η ∂ ∂=−+Φ− ∂Φ ∂ −+ ∂ ∂ +∇ ∗∗ tjj R Rppj , ,, A e e hyp 2 (5) v() () → ρ ρρρ ∂ ∂ =−∇⋅ +∇⋅∇++ ∇ ρ⊥⊥ t DSD , tot hyp 4 (6) vv () ∥∥ →→ γκκ ∂ ∂ =− ⋅∇ −∇⋅+∇⋅ ∇+∇+ ⊥⊥ p t pp TT S, EE T (7) where the density, temperature and current sources ρ S , ST and jA have been introduced. The current source term jA also includes the time-dependent bootstrap current calculated using Sauter’s formula [20]. The convective derivative, the parallel gradient, the perpendicular gradient, and the Poisson brackets are defined as v → = ∂ ∂+⋅ ∇ tt d d , E (8) [] ∥ →→ ∇= ⋅∇bb , (9) ∥ ∇=∇−∇ ⊥ , (10) [] () → αβ αβ=⋅∇×∇ φ e,, (11) →→ = b BB 1. (12) Note that equations(3) and (4) can be reduced to scalar equationsby projecting them in the poloidal and parallel directions, respectively, by applying the operators [ ()] → ∇⋅× φ Re and () →⋅b . This reduced set of equations(without the diffusive transport terms and the diamagnetic terms) is equivalent to that derived in [21], where energy of the system is shown to be conserved at first order. Nucl. Fusion 57 (2017) 076006 S.J.P. Pamela etal 3 The perpendicular mass and thermal diffusivities ⊥ D and κ ⊥ used in simulations are ad hoc coefficients with a well at the pedestal region to represent the H-mode transport barrier. The profiles used for ⊥ D and κ ⊥ are discussed in section2.2. A Spitzer-like resistivity (/ )/ η η=− TT o oee, 32 is used, with To e, the electron temperature at the magnetic axis. Likewise, a temperature-dependent perpendicular viscosity is used: (/ )/ µ µ=− TT ooee,32 . The Braginskii parallel thermal conductivity ∥ κ is expressed as (/) ∥∥ / κ κ=TT oo 52 . The ratio of specific heat is /γ=53 . Hyper-diffusive coefficients µhyp , ηhyp and Dhyp are also used in these simulations. The normalization of the equationsis based on the magnetic permeability µo and the core density ρo , so that time is normalized to a near Alfven time /µρ=t t oo SI . For a deuterium plasma with particle density m=× − n610 o19 3 , a normalized time unit corresponds to approximately 0.5μ s. Naturally, current is normalized with µo and density with ρo . Pressure is also normalized with µo , the diamagnetic frequency is normalised with time, as and δδµρ= ∗∗ oo SI and the diffusive parameters are normalized as /η ηρµ= oo SI , /µ µµρ= oo SI , µρ=DD oo SI and /κκµρ=oo SI . The boundary of the computational domain in the SOL is a flux surface, on which Dirichlet boundary conditions (zero perturbation) are applied for all variables, except density and temperatures, for which Neumann conditions with null gradient are applied. At the divertor targets, sheath boundary conditions are used for the parallel velocity and the energy conduction, such that, γ⋅=±=±bc T , stot v → → (13) vv ∥∥∥∥ →→ κγ+∇=nT TnT. sh (14) where → b is the unit vector parallel to the magnetic field. The density has free outflow boundary conditions at the target (no density reflection). In the private region, which is also bounded by a flux surface, Dirichlet conditions are used for all variables. The 2D poloidal grid is composed with isoparametric cubic Bezier finite elements [4]. The finite element grid is aligned to equilibrium flux surfaces for the three regions of the core, the SOL and the private region. Alignment along open flux surfaces in the SOL is important in order to treat accurately the fast parallel transport of energy along magnetic field lines. The toroidal dimension is represented by a Fourier series. The time stepping is done using the implicit Crank– Nicolson scheme, so that the size of time steps depends only on the time scale of the instabilities that are simulated. This implicit scheme results in a sparse system of equations, which is solved using a Generalized Minimal REsidual Solver (GMRES). The preconditioner for this iterative GMRES is obtained by solving independently each sub-matrix corresponding to different Fourier harmonics, which amounts to a block-Jacobi preconditioner. These sub-matrices are solved using the direct parallel sparse matrix solver PaStiX [22]. In order to allow the n = 0 component of the →→ ×E B and parallel flows to evolve towards a stationary equilibrium, the simulations are first run without toroidal modes, with only the equilibrium n = 0, for 0.5 ms. This allows the Bohm boundary conditions to diffuse into the SOL (at the zeroth time step, v∥ → is Mach-1 on the target, and zero inside the plasma, already at the nodes adjacent to the boundary) 2.2. Simulating ELMs in JET The pulses used for the simulations are the same JET-ILW as in [14]: an IP and a ν* scan, with both lowand highgas fuelling. In the experiments, the divertor heat-flux is observed to increase at lower pedestal collisionality, but one of the most notable effects is that of the gas fuelling: at similar νped * values, a higher gas fuelling leads to higher ELM frequency, with lower ELM sizes and divertor heat-fluxes. In this previous study [14], the ELM energy losses were reproduced by the simulations, but the divertor heat-fluxes were lower than in experiments. The simulated ELMs were less intense, but lasted longer than in experiments (hence ELM energy losses were comparable). In the present study, more accurate pre-ELM equilibria were used, and larger toroidal resolutions as well as diamagnetic terms were included in the simulations. In [14], the standard EFIT equilibria were used, with magnetic constraints only, in order to map the ne and Te profiles and reconstruct the JOREK Grad–Shafranov equilibria. A recent version of EFIT++/JEC2020 [23] was used to impose additional constraints (pressure, polarimetry, motional stark effect MSE) to the pre-ELM equilibria. The most important difference between the two versions of EFIT, with respect to ELM simulations presented here, is that the latter version allows for pedestal currents. This modifies the magnetic flux gradient in the pedestal region, such that when the highresolution Thomson scattering (HRTS) ne and Te profiles are mapped from real space onto the magnetic equilibrium, the Figure 1. The pressure, polarimetry and MSE constraints in EFIT, together with allowing for pedestal current in the equilibrium reconstruction, have a strong effect on the mapping of the ne and Te profiles from the HRTS diagnostic. This plot shows the relative difference, in %, between the standard EFIT and EFIT++, for the pedestal width of the ne, Te and pe profiles. The ψ-mapped pedestal width can be >30% smaller with EFIT++. Nucl. Fusion 57 (2017) 076006 S.J.P. Pamela etal 4 resulting pedestal gradient width (in ψ-space) is modified. For these pulses it is found, as shown in figure1, that the gradient width can be diminished as much as 40% for low collisionality pulses (i.e. those with the largest pedestal current). This effect was expected to have a strong influence on the discrepancy obtained in the first study of [14]. The pre-ELM ne and Te profiles used for the simulations are taken from the high-resolution Thomson scattering (HRTS) diag nostic, as described in [14, 24, 25]. In all simulations, Ti = Te is assumed, which will be discussed further in the context of lowand high-gas fuelling discharges. The pedestal current density is calculated according to the Sauter’s bootstrap cur rent model [20]. For high ν∗ pe d discharges, it has been shown in [26, 27] that the Sauter formula [20] overestimates the bootstrap current, however in the present simulations this discrepancy was not accounted for and the Sauter formula was used for every discharge. Note that the pre-ELM ne, Te and current profiles are not taken directly from the EFIT++ equilibrium itself, because EFIT++ will find the optimal equilibrium that satisfies all various constraints (pressure, polarimetry, MSE and magnetics). This can result in output profiles that are slightly different from the input measured HRTS profiles. Therefore, the original pre-ELM profiles are used, with the consistent bootstrap current, and the Grad–Shafranov equilibrium is recalculated internally in JOREK, as described in [28]. Having equilibrium profiles and currents in coherence with the background ψ map is important in nonlinear simulations, where the ∇p and →→ ×JB terms need to be properly balanced in the momentum equation(3). In simulations below, the full diamagnetic effects are used, as well as multiple toroidal mode numbers, which was not the case in the previous study [14]. The new simulations were run with toroidal mode numbers n = 3,6,9,12,15. A few simulations were run with the full toroidal spectrum from 1 to 15, n = 1, 2, 3, ..., 15 (although this was not possible for all pulses). All simulations were run with a resistivity level at a factor 10 above the Spitzer resistivity, which depends on the absolute temperature profile of each pulse. The viscosity level is kg ms⋅⋅ −−− 4.10811 for all pulses. The parallel conductivity, which also depends on the absolute temperature for each pulse, was taken at a factor 4 above the ion Braginskii value (10x below the electron value). The factor 10 in resistivity results from numerical limitations: at low resistivity, diamagnetic effects are challenging due to strong diamagnetic currents obtained during the ELM crash. To avoid thin current layers below the grid resolution and below the electron gyroradius, hyper-resistivity is set to η=− 10 hyp 13 , which is approximately η×10 2 in the pedestal region, to ensure that MHD is not dominated by the numerical resistivity, rather than resistivity itself. The same value is used for the hyperviscosity µhyp and the hyper-diffusivity Dhyp. The perpendicular diffusivity profiles ⊥ D and κ ⊥ are based on the initial density and temperature profiles in the edge region, with /ρ∼∇ ⊥ D1 and /κ∼∇ ⊥T1 , from ψ=0.5 n up to outside the separatrix (at the bottom of the profiles). Hence, the perpend icular diffusive fluxes ρ∇ ⊥ D and κ∇ ⊥T are constant in the pedestal, which ensures a rigid evolution of the density and temperature profiles. The value of ⊥ D and κ ⊥ is taken to be ms − 521 and ms(  )⋅ −− 10 71 respectively, at ψ=0.8 n , such that the diffusive evolution of the profiles is much longer than the typical time scale of an ELM (∼1–5 ms). Radially uniform sources ρ S and ST are used inside the separatrix to ensure that the (also radially uniform) perpendicular diffusive fluxes are balanced at equilibrium. This particular choice of diffusive profiles and sources was made because it is the most robust way to ensure that the profiles remain identical to the initial profiles, as measured by the HRTS diagnostic. Note that since the simulations are started with only n = 0 (to allow for the formation of the stationary equilibrium), this ensures that the profiles do not deviate from the pre-ELM HRTS profiles in this period, before the n > 0 toroidal modes are started, and the ELM crash begins. Eventually, integrated simulations should be run with more realistic transport and heating/fuelling models to determine what the exact values of diffusive parameters and sources should be in each JOREK simulations. However, such an inclusion of JOREK into integrated modelling structures would require a significant investment which, at this point, cannot be characterised as essential with respect to the peeling–ballooning instabilities themselves. In any case, there is no robust and reliable models for turbulent transport in the pedestal region in the presence of an H-mode transport barrier. The grid chosen for the simulations has elements radially concentrated in the pedestal region, so that the radial width of elements is 1.4 mm at the separatrix (at the midplane). The radial resolution in the core is kept high enough to avoid artificial ballooning instabilities with the highest mode number n = 15. The poloidal distribution of elements is also controlled to ensure a uniform poloidal resolution, including at the X-point. The average poloidal width is 2.7 cm . Since the time step of implicit schemes depends only on the nonlinear activity of the toroidal modes (i.e. not on the minimum element size of the 2D grid), the usual time step of the order of 0.5 μ s , which is reduced up to 0.05 μs in the most challenging cases. This means that a usual ELM simulation requires between 2000 and 10 000 time steps. 3. Simulation results for ILW pulses 3.1. Linear stability The first aspect of simulations to consider is the linear stability of the simulated discharges. It has now been established that, in JET-ILW type-I ELMy discharges with high gas fuelling, the linear stability of pre-ELM equilibria is typically incoherent with the linear ideal MHD stability theory, unless Nitrogen seeding is used [26, 29]. Linear MHD calculations of pre-ELM stability with the ELITE code [30, 31] shows that the ELM onset is typically well inside the stable region of the j-α stability diagram, where j is the edge plasma density (responsible for destabilizing peeling modes), and α is the pedestal pressure gradient (responsible for destabilizing ballooning modes). The linear stability of peeling–ballooning modes can also be evaluated using JOREK, including non-ideal MHD effects, such as resistivity, viscosity, diffusion and conductivity, NBI Nucl. Fusion 57 (2017) 076006 S.J.P. Pamela etal 5 toroidal rotation, and diamagnetic effects. This is achieved by starting simulations with a stable pedestal pressure, progressively increasing it through the stability threshold, as shown in figure2(a). Once a coherent peeling–ballooning mode starts growing exponentially, the pedestal pressure level is used as the stability threshold. Note that this is done using a fixed pedestal width, and with density and temperature sources that result in the same ne ped and Te ped values as in the experiments once the experimental pped is reached. Also, the bootstrap current (calculated using Sauter’s formula [20]), is evolving according with the pedestal density and temperature profiles. The results for linear ideal and non-ideal MHD calculations are shown in figure2. Figure2(a) shows how JOREK is used to evaluate the linear stability of the pedestal by progressively increasing pped at fixed pedestal width. Figure2(b) shows an example of the j-α diagram for the JET-ILW pulse 83340, including the experimental pre-ELM equilibrium, the ideal MHD threshold (with the stability boundary) calculated by ELITE, and the non-ideal MHD threshold calculated by JOREK. Figure2(c) shows the ideal and non-ideal calculations for all the ILW pulses considered in this study. This shows that in many cases the non-ideal MHD effects have a significant impact on the linear stability threshold. At this point it is not clear which of these particular effects has the most important impact on stability, and further work would be required to isolate each effect. It will be the scope of a future study, since this requires re-running each discharge 5 to 6 times (once for each isolated effect), which is not achievable on a short time scale. This result demonstrates that JOREK could be used to improve our undestanding of MHD stability in JET-ILW (a) (b) (c) Figure 2. (a) A JOREK simulation starting from a stable pped value. The top plot shows the evolution of pped, the bottom plot shows the evolution (on a logarithmic y-axis) of the kinetic energy of various ballooning modes. Keeping the pedestal width constant and increasing pped progressively reveals the linear MHD threshold for peeling–ballooning modes (determined once a mode starts growing exponentially, here at ∼t1800 ). The advantage over linear MHD codes is that non-ideal effects (resistivity, rotation etc) are easily included. The disadvantage is that it is numerically more expensive. (b) An example is shown here, for pulse 83340, to illustrate how the JOREK code can, in some cases, obtain a better agreement with experiments than ideal MHD calculations. This is shown in the usual j-α diagram (the solid red line shows the stability boundary calculated by ELITE). The dashed blue line represents the coherent evolution of αmax and jmax, assuming that the bootstrap current evolves when pped is increased (here according to the Sauter model). (c) This plot shows the linear calculations for all JET pulses included in this study, with the ELITE code (for ideal MHD), and with the JOREK code (for nonideal MHD). These calculations are compared to the experiments on the x-axis. Non-ideal MHD clearly provides a closer agreement with experiments. The red square/star is pulse 83340. Nucl. Fusion 57 (2017) 076006 S.J.P. Pamela etal 6 discharges. In particular, once the physics responsible for this correction on ideal MHD is isolated, JOREK could be used to provide a reinforced MHD stability calculation for peeling–ballooning modes in predictive models like EPED [32]. However, at this point a systematic inclusion of JOREK calculations into predictive models is strongly limited by the fact that such JOREK calculations (as shown in figure2(a)) are numerically more expensive than ideal MHD calculations using dedicated codes like ELITE. 3.2. ELM Size and divertor heat-fluxes As shown in figure2, linear non-ideal MHD provides a better agreement with experiments than ideal-MHD for linear calcul ations. However, linear stability provides no estimate or guarantee concerning the size of the ELM crash itself. As contradictory as it may seem, while reasonable agreement is obtained for linear stability, the ELM crash itself is typically smaller in simulations. However, it is worth noting that the discharges which give the best agreement (in terms of ELM size and divertor heat-flux) are those closest to the linear ideal MHD threshold for PB-stability, which are typically the low gas fuelling discharges. In this section, the ELM simulations are obtained by starting from the exact pre-ELM equilibrium prescribed by experimental data. (i.e. not from a progressively increasing pressure gradient as described above.) The ELM size calculated in JOREK simulations with diamagnetic terms are smaller than in the experiments by a factor 2–3. This is a large decrease compared to the previous set of simulations without diamagnetic effects, as shown in figure3(a). The experimental ELM size is determined using the preand post-ELM pressure profiles obtained from the HRTS diagnostic, as in [14]. It should be noted that the diamagnetic effects are dominating the non-linear phase of the ELM crash, rather than the linear burst and the initial filamentation of the plasma. Of course, linear stability is also affected, to the extent that some pulses that are PB-unstable without diamagnetic terms become stable with diamagnetic terms. In most cases however, with the diamagnetic terms, the initial burst of filaments is similar, but the duration of the instability is much shorter with diamagnetic terms. Without ω* effects, filaments keep bursting through the separatrix until the excess pressure is evacuated from the pedestal. However with ω* effects, at most one or two sets of filaments burst through the separatrix, after which the stabilising effect of the ω* terms set an end to the ELM crash. This aspect of the shortened ELM duration is certified by the dynamics of divertor heat fluxes. Figure 3(b) shows the peak heat-flux on the outer divertor for each pulse, as a function of pedestal collisionality ν∗ pe d . The experimental data is obtained using the Infra-Red camera diagnostic [1, 33], averaged over all type-I ELMs for each discharge. The divertor heat-flux is similar, for most pulses, to the previous simulations of [14]. The width of the outer divertor heat-flux (10–12 cm) is also similar to the experiments and to the previous simulations. The main effect of diamagnetic terms is to diminish the time duration of the ELM crashes. Note that at the lowest ν∗ ped , there is a large span in peak divertor heat flux for the experimental points due to different gas fuelling levels, which has a strong influence on the ELM frequency, and therefore the ELM amplitude. It should be pointed out that, of course, although simulations without diamagnetic effects give a better agreement with the experimentally measured ELM size, they are not more accurate than simulations without diamagnetic terms. On the contrary, the increased duration of the ballooning activity without diamagnetic terms should be seen as a physical artefact, even if it leads to the right ELM energy losses. In fact, the simulations without diamagnetic terms are more coherent in the sense that both the divertor heat-fluxes and the ELM size are smaller than in the experiments. Therefore, the issue is Figure 3. (a) This plot shows the total ELM energy losses as a function of pped, for the experiments (red circles, calculated using preand post-ELM HRTS profiles), for the simulations without diamagnetic terms (black stars), and for simulations including diamagnetic terms (blue triangles). The diamagnetic effects stabilise the mode activity in the late phase of the ELM, leading to shorter bursts, and hence smaller ELM energy losses. (b) The peak heat flux on the outer divertor as a function of ν∗ pe d , for the experiments (red circles), for the simulations without diamagnetic terms (black stars), for simulations including diamagnetic terms (blue circles), and for two theoretical ν∗ pe d scans based on pulses 83330 and 83334 (green triangles). Nucl. Fusion 57 (2017) 076006 S.J.P. Pamela etal 7 now to determine why the divertor heat-flux is smaller, since it represents the signature of up-stream MHD activity. To address this question, we consider the cases that do have a good agreement with the experiments. In particular, two cases are relatively well reproduced, both in terms of divertor heat flux and energy losses: discharges 83330 and 83334. These are low gas fuelling discharges with MA TMW[][ ]=IBP,, 2,2,20 PTNBI and MATMW[][]=IBP,, 2.4,2.4,25 PTNBI respectively. Using these two cases, a theoretical ν∗ pe d scan was performed by varying the density and temperature levels at constant pressure. This ensures that the ideal MHD ballooning stability remains identical. In practice, this is equivalent to changing the normalisation of the non-ideal MHD parameters, such as resistivity, viscosity, thermal conductivity, and diamagnetic effects (as described in section2.1). The result is in good agreement with the experiments, as shown in figure3(b). This is important as it demonstrates that the parallel transport of energy along field lines is relatively well described by the Braginskii model. Another quantity to consider for the validity of parallel transport with a Braginskii model is the parallel energy, as defined by Eich etal [34]: MJ m()[] () ∥∥ ∫ ε α⋅= − sq st t,d t B 2 ELM (15) where s is the divertor (radial) coordinate, and αB is the angle between magnetic field lines and the divertor target. The regression scaling obtained from the experimental data on multiple tokamak devices [34] shows a strong dependency on pped and on the relative ELM size ∆WELM : ε ⋅=|× × ×∆ × −± ±± ±± nT WR 0.28 . 2 0.14 e,ped 0.75 0.15 e,ped 0.98 0.1 ELM 0.52 0.16 geo 1.00.4 MJ m[] ∥ (16) Both (15) and (16) can be calculated in simulations, and it is found that for cases with an ELM size that corresponds to the experiments, the Eich scaling is well reproduced. This is shown in figure 4, where the two formulas are compared, both for the cases with and without diamagnetic terms. The agreement is successful for the cases without diamagnetic terms, while there is a larger disparity for the cases with diamagn etic terms, due to the discrepancy in ELM size. However, the parallel heat transport model is the same regardless of the diamagnetic terms, and therefore, if future simulations with diamagnetic terms can be improved to match the exper imental ELM size, this indicates that the Eich scaling should be recovered. Hence, if the linear stability from the experiments is well reproduced, and the parallel energy transport is also coherent with experimental regressions (provided the ELM size is reproduced), then the next logical step would be to look at the physics model. In particular, since low-gas pulses like 83330 and 83334 are relatively well reproduced by JOREK, the neutrals model [35, 36] could be considered in future simulations. However, it should be noted that one effect from neutrals can already be eliminated, namely the outward density shift observed at high gas fuelling [37–39], since this is already taken into account in the pre-ELM equilibria. Hence, if the missing neutrals were responsible for the lower ELM energy losses in simulations, it would have to do with other effects. For example, at high gas fuelling, the temperature conducted along field lines would face larger parallel gradients in the divertor region due to stronger ioniz ation and charge exchange near the target, which would increase the parallel energy fluxes ∥∥ κ∇T in the divertor region and thus increase the rate at which the pedestal temper ature is evacuated. Another example is that the shift observed between the pre-ELM ne and Te pedestal profiles at high fuelling should be verified for the ion temper ature Ti as well, which is not considered in present simulations, due to the time-resolution of the charge exchange diagnostic for these pulses (i.e. the Ti profiles are ELM-averaged). Hence, if the Ti profile shifts together with ne (as opposed to having =TT ie ), the stability of the highfuelling discharges would be strongly affected. Such questions should be addressed in future studies. Figure 4. (a) The parallel energy arriving on the outer divertor can be compared against the Eich scaling law equation(16). The agreement is good provided the ELM size is comparable to the experiments. Hence, the best agreement is for the cases without diamagnetic terms. Note that the Eich scaling used here is the exact regression ε=∆ nT WR 0.28 e,ped 0.75 e,ped 0.98 ELM 0.52 ge o 1.0 ∥ . (b) The same plot, but using modified exponents, within the regression limits: ε=∆ nT WR 0.14 e,ped 0.6 e,ped 0.98 ELM 0.68 ge o 1.0 ∥ . Nucl. Fusion 57 (2017) 076006 S.J.P. Pamela etal 8 3.3. Nonlinear stability and multi-ELM cycles There is, however, another aspect of simulations which should be addressed, because it could have a significant impact on ELM energy losses: the nonlinear dynamics of ELMs, and more particularly, the nonlinear stability of ELMs. One of the peculiar aspects of the simulations presented in section3.2 is the absence of nonlinear coupling between modes during the ELM crash. In most cases, a dominant mode number leads to the crash, while minor coupling with the other modes occurs only at the end of the ELM, once most of the energy has already been evacuated. In few cases, the coupling occurs earlier in the crash, but there is usually still one dominant mode number. It has been experimentally established that large type-I ELMs have a strong nonlinear comp onent [40–42], notably for the main part of the crash. Thus, the quasi-linear aspect of ELMs in the simulations, even in the presence of diamagnetic terms, suggests that the set-up of simulations should be reconsidered. In practice, any initial-value code like JOREK needs to start the non-axisymmetric perturbations as low numerical noise, which eventually organises into coherent peeling–ballooning structures, provided the equilibrium is unstable with respect to these instabilities. These peeling–ballooning structures will then grow exponentially, until they reach an amplitude large enough to perturb the background axisymmetric equilibrium, and thus create a crash. However, each toroidal Fourier harmonic, which represents an individual mode number, will grow at its individual rate. Since these growth rates depend on the mode number, one mode will grow more rapidly, and hence reach the equilibrium before all other modes. This is shown in figure5, which was run with the full spectrum of modes n = 1, 2, 3, ..., 15. This simulation was performed to test whether one-to-one coupling between modes was required to obtain more nonlinear interactions between modes in the early phase of the ELM. However, the result is similar to the filtered spectrum n = 3, 6, 9, 12, 15: nonlinear coupling only occurs in the late phase of the ELM. It is worth noting that the level of coupling in the early phase of the ELM strongly depends on the choice of MHD parameters, particularly resistivity and viscosity. In cases with higher resistivity and viscosity, nonlinear coupling in the early phase of the ELM is frequently observed, as in [17, 18]. However, the theoretical resistivity is set by the Spitzer value, and a higher viscous coefficient would stabilise the ballooning modes, which is the reason for the present choice of η and μ. In experiments, the whole spectrum of modes is active at a certain level, not far below the equilibrium, and thus when an ELM arrives, all modes will crash together. It could even be assumed that nonlinear coupling is occurring before the ELM crash itself. In fact, it seems that the nonlinear coupling itself is responsible for the crash, as proposed in [43]. In order to simulate this, a state of marginal instability must be achieved, in which the various toroidal modes are fluctuating at equilibrium level without causing a crash. This is equivalent to going progressively from a stable to an unstable equilibrium, through the ELM onset threshold. In other words, it is equivalent to simulating multiple type-I ELM cycles. JET pulse 83334 is started from a stable pedestal pressure, which is then restored using density fuelling and heating, a first small quasilinear crash is obtained, similar to those presented in figure 3. Further increasing the pedestal pressure beyond this initial crash, strong nonlinear coupling is obtained, leading to a second large pedestal crash. This is represented in figure6(a), which shows the evolution of the (a) (b) Figure 5. (a) The kinetic energy of the modes n = 1, ...,15 (in logarithmic scale) as a function of time. The mode n = 15 has a slightly higher growth rate than the other modes. However, since the perturbations are started at a low numerical noise level, the exponential growth phase is long enough to cause a large delay between the arrival of n = 15 and n = 14 at equilibrium level. Therefore, by the time n = 14 reaches equilirbium, n = 15 has already produced a crash, which reduces pped, stabilising all other modes, and thus preventing any nonlinear coupling between modes. Note: the discontinuities observed at 0.40 and 0.43 ms are just a jump in numerical noise (those toroidal harmonics have not yet organised into peeling–ballooning structures). (b) The same plot, but not in logarithmic scale, and zoomed on the ELM crash itself, to show the obvious absence of any nonlinear coupling. Nucl. Fusion 57 (2017) 076006