Sylvain Bertho

Radhouan Maazaoui

Damla Torun

Ismael Dondasse

Raoudha Abderrahim

Isabelle Gillaizeau

Cyril Nicolas

Ismaël Dondasse

Iron Catalyzed B-C(sp2 )-H Alkylation of Enamides

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.

Introduction

Enamides are highly valuable structural motifs that are gaining increasing interest, due to their diverse possible transformations, and also because they are found in many natural products and pharmaceutical drugs (Figure 1). [1][2][3][4][5][6][7][8][9][10] As a result, over the past decades, various synthetic methods have been considered on enamide substrates. 11 In recent years, a special focus has been placed on the development of their direct -C(sp 2 )-H functionalization, including arylation, olefination, alkynylation, carbonylation, sulfonylation and phosphorylation reactions among a few other useful transformations. 11b,11e-11i Remarkably, alkylation and/or fluoroalkylation reactions, which represent a greater challenge as they involve a β-C(sp 2 )-C(sp 3 ) bond formation, have also been at the forefront of investigations. 11b,11f,11i Scheme 1. Strategies for direct intermolecular (fluoro)alkylation of enamides.

Significant advances have been achieved in this area, mostly through the addition of alkyl radicals, which are mainly generated by direct radical initiation or a single electrontransfer (SET) process using copper catalysts or iridium(III) salts under photoredox-catalysis (Scheme 1a). 11b From an environmental and economical point of view, iron catalysis holds the promise of being able to cover a wide swath of organic synthesis. 12 Iron catalysts are not only expected to operate through reductive and oxidative regimes, but they are also able to embrace "early" transition metal characters. Of note, in 2013, Guan and co-workers developed a FeCl3-catalyzed selfcondensation of enamides, offering an original strategy for the synthesis of nitrogen-containing quaternary carbon centers through Lewis acid catalysis (Scheme 1b). [START_REF] Zhao | [END_REF] In 2015, our group reported the first example of the iron(II)-catalyzed direct C-H trifluoromethylation of enamides using Togni's reagent via a radical process (Scheme 1c). 14 Fluorine is indeed an increasingly commonly found component in bioactive compounds that span all therapeutic categories. Fluorine-containing derivatives are amongst some of the top-selling and best-performing "small molecule" pharmaceutical products. [START_REF] Hiyama | Organofluorine Compounds: Chemistry and Applications[END_REF][START_REF] Inoue | [END_REF] However, to the best of our knowledge, there are no reports so far on the use of iron catalysis with the perspective of developing an efficient, straightforward and mild application with broad reagent scope in order to functionalize enamides through a direct intermolecular -C(sp 2 )-H (fluoro)alkylation. Thus, there is still great room for improvement. We report herein our efforts to overcome this challenge.

Results and Discussion

In our endeavors to explore the reactivity of iron catalysis in the alkylation reaction of enamides, we started our study with the six-membered cyclic enamide 1a and the ethyl 2-bromo-2,2difluoro-acetate 2a (2 equiv) which were selected as model substrates (Table 1). Based on the conditions previously developed by our team for the copper-catalyzed olefinic C-H difluoroacetylation, 17 and iron-catalyzed trifluoromethylation of electron-rich olefins, 14 the optimization of the reaction conditions was started by using iron(II) dichloride (0.2 equiv.) as the catalyst, 1,10-phenanthroline (0.24 equiv.) as the ligand, and K2CO3 (2 equiv.) as the base. The reaction proceeded at 100 °C in acetonitrile for 16 h under an argon atmosphere. The desired C3 alkylated enamide 3aa was obtained as the unique derivative, with complete regioselectivity in 42% yield (entry 1). To improve the yield of 3aa, screening of the catalyst loading (entries 2-3) and of the reaction time (entry 4) was thereafter carried out with success with yields up to 91%. The efficiency of the amount of α-bromo reagent (entry 5), base (entries 6-8), and solvent types (entries 10-14), as well as iron (entries 15-21) and ligand source (entries 23-24) were also evaluated, but no improvement whatsoever of the reaction outcome was observed. Importantly, when no base (entry 9) or no catalyst (entry 22) was used, the alkylation was not promoted, and no product was noted. Iron complexes of η-nitrile ligands such as acetonitrile are well known systems capable of mediating various useful synthetic transformations. [START_REF] Stephenson | Knowledge Updates 2014[END_REF] Of note, when the alkylation reaction was performed in acetonitrile, in the absence of 1,10-phenantroline, the enamide product was provided, albeit with lower efficiency (31% vs. 86%; entry 25 vs. entry 5). Interestingly, when iron(II) or iron(III) salts were employed, the reaction was also amenable to compound 3aa (entry 21). a Reaction conditions: 1a (0.5 mmol, 1.0 equiv.), 2a (2-3 equiv.), 1,10phenanthroline (0.12 mmol, 0.24 equiv.), dry base (1.0 mmol, 2.0 equiv.), solvent (3.0 mL), 100 °C, 16 h, under Ar. b Yields determined by [START_REF] Shen | [END_REF] F NMR analysis of the crude mixture using α,α,α-trifluorotoluene as an internal standard. To our delight, most of the reactions proceeded smoothly to afford the corresponding -alkylated enamides 3 in moderate to good yields. A good functional group tolerance was observed except for the substrates 1b and 1d which are more electronrich than the other enamides because of the presence of the N-PMB protective group (1b) or of the oxygen atom in the morpholine moiety 3d. Interestingly, in the difluoro series, the deactivated N-CO2Ph-protected piperidinone adduct 3ca was isolated in moderate yield as a 6:4 ratio of rotamers. Of note, a C-5 selective alkylation of 2-pyridone 1e and 1f was observed with good yields, and no trace of the other regioisomers was detected. Reaction with N-acyl enamide 1i was also possible and a decent 49% isolated yield of 2ib was achieved. In addition, the acyclic enamides 1g-h were well tolerated and a clean regio-and excellent diastereoselectivity was observed (e.g., 3gh). Their E configuration was confirmed by nOe and 1 H-NMR experiments (see SI). With a view to enhancing molecular diversity, the reagent scope was exemplified onto exo-enamides 1g-h and a few other derivatives in the piperidinone and morpholinone series. Remarkably, similar results were obtained for all derivatives, leading to the corresponding alkylated adducts in moderate to good yields (59-83%), apart from compounds 3ae, 3bc and 3dc which were respectively obtained in non-optimized 44%, 42% (63% brsm) and 25% (47% brsm) isolated yields.

Scheme 2. Mechanistic studies on FeCl2-catalyzed alkylation of enamides 1.

To have a better understanding of the Fe(II)-catalyzed -C(sp 2 )-H alkylation of enamides 1, some mechanistic experiments were carried out. As depicted in scheme 2 (a), when the reaction of 1a and 2a was conducted under the above optimal conditions in the presence of 3 equiv. of TEMPO (2,2,6,6-tetramethyl-1piperidinyloxy) free radical, the reaction was suppressed and the desired product 3a was obtained only in trace quantity. A complete inhibition was also observed in presence of BHT (2,6di-tert-butyl-4-methylphenol). These results indicate that the reaction might proceed through a radical pathway. Based on these experimental results and previous reports, 14,[START_REF] Shen | [END_REF] Please do not adjust margins Please do not adjust margins valent iron complex, Fe II is regenerated leading to the formation of a carbocation species C and its corresponding resonance structure C'. Finally, a base mediated deprotonation furnishes the desired alkylated product derivative 3 with a total transdiastereoselectivity (R 1 ≠ O).

Finally, this chemistry can be advantageously extended to lowcost and high-abundant iron powder. Noteworthy, we noticed a good conversion of the acyclic enamides of type 1g and 1h when treated with difluoroacetate 2a and a stoichiometric equivalent of iron-powder in DMF at 110 °C. Interesting results using iron powder were recently described by Iwasaki, Nishihara and co-workers on a practical synthesis of γ-lactones from alkenes and α-halocarboxylic acids and their derivatives. 20 As depicted in Scheme 3, iron powder showed the catalytic activity needed to drive the C(sp 2 )-H alkylation of enamides 1g-1h in the presence of α-haloesters 2a-d (e.g., G.P. B). To our delight the reaction proceeded smoothly without using ligands, bases or additives to afford the related adducts in good yields, with again a total control of the regio-and stereoselectivity in favor of the E-diastereomers. However, surprisingly, this alkylative coupling was not amenable to cyclic endo-enamides.

Conclusions

In summary, we have developed a mild and novel method for the -C(sp 2 )-H (fluoro)alkylation of cyclic and acyclic enamides using either a catalytic or a stoichiometric amount of iron species and various α-bromo activated reagents. Mechanistic studies reveal that the reaction might proceed through a radical process. Interestingly, the -C(sp 2 )-H alkylation of enamides was also originally achieved using a stoichiometric amount of low-cost and high-abundant iron powder without using additional ligands, bases, or additives. Further applications based on this chemistry are in progress in our laboratory.

Experimental section

singlet, d = doublet, t = triplet, q = quartet, m = multiplet, o = overlapped and po = partially overlapped), coupling constant Ј (Hz), and integration. [START_REF] Zhao | [END_REF] C NMR (63, 101 or 176 MHz) chemical shifts are given in ppm. Spectra were calibrated using the corresponding nondeuterated solvent. [START_REF] Shen | [END_REF] F NMR (235 or 376 MHz) chemical shifts are given in ppm. High-resolution mass spectra were recorded with a Bruker maXis ESI qTOF ultrahigh-resolution mass spectrometer coupled to a Dionex Ultimate 3000 RSLC system. MS data were acquired in positive mode and were processed using Data Analysis 4.4 software (Bruker). Infrared spectra were recorded neat with a Thermo Scientific Nicolet IS10 FTIR spectrometer using diamond ATR golden gate sampling and are reported in wave numbers (cm -1 ). Analytical thin-layer chromatography (TLC) was performed with Merck Silica Gel 60 F254 precoated plates. Visualization of the developed chromatogram was performed under ultraviolet light (254 nm) and on staining by immersion in aqueous, acidic ceric ammonium molybdate followed by charring at ca. 150 °C. Column chromatography was performed in air on Silica Gel 60 (230-400 mesh) with petroleum ether (PE, bp 40-65 °C) and ethyl acetate as eluents, unless otherwise stated. Organic solutions were concentrated under reduced pressure with a Buchi rotary evaporator. The IUPAC name of the new compounds was generated automatically using the structure-to-name generator included in BIOVIA Draw 2020.

General Procedure for the β-C(sp 2 )-H Alkylation of Enamides 1a-i Using Catalytic FeCl2 (G.P. A). An oven-dried microwave vial under argon atmosphere was charged with the enamide substrate 1 (0.5 mmol, 1.0 eq.), FeCl2 (13 mg, 0.1 mmol, 0.2 eq.), 1,10-phenanthroline (22 mg, 0.12 mmol, 0.24 eq.), dry K2CO3 (138 mg, 1.0 mmol, 2.0 eq.) and a magnetic stir bar. Anhydrous acetonitrile (3.0 mL) was added and the flask was submitted to three freeze-pump-thaw-argon cycles and placed under argon atmosphere. Next, the corresponding αbromo reagent 2 (1.5 mmol, 3.0 eq.) was added, the reaction vessel was sealed and the reaction mixture was heated at 100 °C for 16 h. After cooling to rt (ca. 20 °C), EtOAc and H2O were then added and

Figure 1 .

 1 Figure 1. Selected biologically active molecules containing enamide functionalities.

of enamides. a, b a

 b Reaction conditions: Enamides 1a-i (0.5 mmol), α-bromo reagents 2a-f (1.5 mmol), 1,10-phenanthroline (0.12 mmol), FeCl2 (0.1 mmol), dry K2CO3 (1.0 mmol), CH3CN (3.0 mL), 100 °C, 16 h, under Ar atmosphere. b Isolated yields. c brsm = based on recovered starting material.

aScheme 3 .

 3 Scheme 3. Substrate scope for the iron-powder-catalyzed -C(sp 2 )-H alkylation of enamides 1ga-hd. a,b

Table 1 .

 1 Optimization studies of the iron(II)-catalyzed -C(sp2)-H alkylation of enamide 1a a

	Entry	Fe source	BrCF2CO2Et	Base	Solvent	Yields
			(equiv.)			(%) b
	1 c	FeCl2	2	K2CO3	CH3CN	42
	2 d	FeCl2	2	K2CO3	CH3CN	79
	3 e	FeCl2	2	K2CO3	CH3CN	65
	4 f	FeCl2	2	K2CO3	CH3CN	91
	5	FeCl2	3	K2CO3	CH3CN	86(82 g )
	6	FeCl2	3	Na2CO3	CH3CN	80
	7	FeCl2	3	Cs2CO3	CH3CN	41
	8	FeCl2	3	Et3N	CH3CN	9
	9	FeCl2	3	-	CH3CN	0
	10	FeCl2	3	K2CO3	CH2Cl2	63
	11	FeCl2	3	K2CO3	DCE h	10
	12	FeCl2	3	K2CO3	Toluene	31
	13	FeCl2	3	K2CO3	DMSO	N.D. i
	14	FeCl2	3	K2CO3	Dioxane	47
	15	FeBr2	3	K2CO3	CH3CN	57
	16	Fe(OTf)2	3	K2CO3	CH3CN	67
	17	FeCl2.4H2O	3	K2CO3	CH3CN	52
	18	Ferrocene	3	K2CO3	CH3CN	19
	19	Fe(OAc)2	3	K2CO3	CH3CN	35
	20	Fe(acac)2	3	K2CO3	CH3CN	3
	21	FeCl3	3	K2CO3	CH3CN	60
	22	-	3	K2CO3	CH3CN	0
	23 j	FeCl2	3	K2CO3	CH2Cl2	21
	24 k	FeCl2	3	K2CO3	CH2Cl2	12
	25 l	FeCl2	3	K2CO3	CH3CN	31
	26 m	FeCl2	3	K2CO3	CH3CN	18
	27 n	FeCl2	3	K2CO3	CH3CN	18
	28 o	FeCl2	3	K2CO3	CH3CN	54

  Please do not adjust marginsPlease do not adjust margins with a range of enamides (1a-i) in presence of α-bromo activated reagents (2a-f). The results are shown in Table2.

	One may therefore wonder about the mechanism and oxidation
	state of the iron catalyst. In fact, it could suggest possible
	Fe II /Fe III catalytic cycle. Lastly, changing the heating
	temperature and the mode, or performing the reaction in open
	air (entries 26-28), led to a poorer performance in the
	transformation, revealing that the best conditions are those
	reported in entry 5.

c 0.1 equiv. of FeCl2 was used. d 0.2 equiv. of FeCl2 was used. e 0.3 equiv. of FeCl2 was used. f Reaction time of 42h. g Isolated yield. h DCE = 1,2-dichloroethane was used. I N.D. = not determined. j 1,2-Bis(diphenylphosphino)ethane, instead of 1,10phenantroline was used as ligand. k N,N,N′,N′-tetramethylethylenediamine (TMEDA), instead of 1,10-phenantroline was used as ligand. l No ligand. m Reaction mixture heated at 80 °C. n The reaction mixture was heated at 150 °C for 1h under μW irradiation. o Reaction performed in open air. With the optimal conditions in hand (i.e., G.P. A), we next explored the scope of this iron-catalyzed C(sp 2 )-H alkylation This journal is © The Royal Society of Chemistry 20xx J. Name., 2013, 00, 1-3 | 3

Table 2 .

 2 Substrate and reagent scope for the FeCl2-catalyzed -C(sp 2 )-H alkylation

  When appropriate or in the event of ambiguous proton and carbon, assignments were established using[START_REF] Shen | [END_REF] F NMR and Heteronuclear Multiple-bond Correlation (HMBC). Unless otherwise stated,[START_REF] Zhao | [END_REF] C and[START_REF] Shen | [END_REF] F spectra were acquired on a broad band decoupled mode.1 H NMR (250, 400 or 700 MHz) chemical shift values are listed in parts per million (ppm), downfield from TMS as the internal standard or relative to the corresponding nondeuterated solvent. Data are reported as follows: chemical shift (ppm on the δ scale), multiplicity

	General Remarks. Unless otherwise stated, all reagents and
	starting materials were purchased from commercial sources and
	used as received. Iron (≥99%, reduced, powder (fine)) and Iron(II)
	chloride (99.99% trace metal basis) were used as catalysts. Toluene
	(ACS reagent, ≥ 99.5%), Acetonitrile (ACS reagent, ≥ 99.5%), and N,N-
	Dimethylformamide (ACS reagent, ≥ 99.8%) were purified by passage
	through a column containing activated alumina under nitrogen
	pressure (Dry Solvent Station GT S100, Glass Technology, Geneva,
	CH). They were subsequently used as solvent in reactions under
	argon atmosphere. Potassium carbonate was activated at 100 °C in
	an oven for 24 h prior to use. NMR spectra were recorded at 298 K
	with a Bruker DPX 250 MHz or Bruker Avance III HD nanobay 400 and
	700 MHz spectrometers equipped with BBO probes. The structures
	of the new compounds were assigned with the aid of 1 D [ 1 H NMR,
	13 C NMR, Distortionless Enhancement by Polarization Transfer
	(DEPT)] and 2 D Correlation Spectroscopy [(1H-1H COSY, and 1 H-13 C
	Heteronuclear Single Quantum Coherence (HSQC)] experiments.

Acknowledgements

The authors thank the CNRS, LABEX SynOrg (ANR-11-LABX-0029), the ANR-20-CE07-0016-01 and the Region Centre-Val de Loire for financial support.

This journal is © The Royal Society of Chemistry 20xx J. Name., 2013, 00, 1-3 | 5

Please do not adjust margins Please do not adjust margins the aqueous phase was extracted with EtOAc (2 ×). The combined organic layers were washed with brine, dried over MgSO4 and concentrated in vacuo. The crude product was purified through silica gel column chromatography to give the desired alkylated enamide product 3 in low to good yield.

General Procedure for the β-C(sp 2 )-H Alkylation of Enamides 1gahd Using Iron Powder (G.P. B

). An oven-dried microwave vial under argon atmosphere was charged with the enamide substrate 1 (0.5 mmol, 1.0 equiv.), Iron powder (58 mg, 1.0 mmol, 2.0 eq.), and a magnetic stir bar. Anhydrous N,N-Dimethylformamide (DMF, 3.0 mL) was added and the flask was submitted to three freeze-pump-thawargon cycles and placed under argon atmosphere. Next, the corresponding α-bromoester 2 (1.5 mmol, 3.0 eq.) was added, the reaction vessel was sealed and the reaction mixture was heated at 110 °C for 4 h. After cooling to rt (ca. 20 °C), the reaction mixture was filtered over celite®, the cake was rinsed (EtOAc) and the solvents were evaporated under reduced pressure. DMF was removed by coevaporating the crude mixture with heptane (3 ×) and the crude residue was purified by flash silica gel column chromatography to give the desired alkylated enamide product 3 in moderate to good yield.

Conflicts of interest

There are no conflicts to declare.


E Nardon
email: eric.nardon@cea.fr

D Hu

F J Artola

D Bonfiglio

M Hoelzl

A Boboc

P Carvalho

S Gerasimov

G Huijsmans

V Mitterauer

N Schwarz

H Sun

Thermal quench and current profile relaxation dynamics in massive-material-injection-triggered tokamak disruptions

3D non-linear magnetohydrodynamic simulations of a disruption triggered by a massive injection of argon gas in JET are performed with the JOREK code. The key role of the thermal drive of the m = 2, n = 1 tearing mode (i.e. the drive from helical cooling inside the island) in the disruption process is highlighted by varying the amplitude and position of the argon source across simulations, and also during a simulation. In cases where this drive persists in spite of the development of magnetic stochasticity, which is favoured by moving the argon source in an ad hoc way from the plasma edge into the 2/1 island at some point in the simulation, a relaxation in the region q ≤ 2 (roughly) takes place. This relaxation generates a plasma current spike comparable to the experimental one. Simulations are compared in detail to measurements via synthetic diagnostics, validating the model to some degree.

Introduction

The possibility of disruptions is a well-identified problem of tokamaks and a threat to the success of ITER [START_REF] Sun | on Disruptions, Plasma Control, and MHD and ITER Physics Basis Editors[END_REF] [2]. In response to this threat, a Disruption Mitigation System ‡ See the author list of M. Hoelzl et al., 'The JOREK non-linear extended MHD code and applications to large-scale instabilities and their control in magnetically confined fusion plasmas', submitted to Nucl. Fusion, preprint at https://arxiv.org/abs/2011.09120 § See the author list of E. Joffrin et al. 2019 Nucl. Fusion 59 112021

(DMS) based on Shattered Pellet Injection (SPI) is being designed for ITER [3][4] [START_REF] Lehnen | 27th IAEA Fusion Energy Conference[END_REF].

In support of the ITER DMS design and in relation with SPI experiments in present tokamaks, SPI-triggered (mitigated) disruptions are being investigated via 3D MHD simulations, in particular with the JOREK [START_REF] Hu | [END_REF][7] [8][9] [10], NIMROD [11] and M3D-C1 [START_REF] Lyons | Verification and Validation of Extended-Magnetohydrodynamic Modeling of Disruption Mitigation[END_REF] codes. Prior to the development of SPI, the main Massive Material Injection (MMI) technique was Massive Gas Injection (MGI), beginning with early experiments at DIII-D [START_REF] Taylor | [END_REF]. MGI was also investigated via 3D MHD simulations, mainly with NIMROD [14][15] [16] but also with JOREK [17][18] and M3D-C1 [19]. Even though the above-cited 3D MHD modelling works validate the codes to some extent, complete quantitative validation has not been demonstrated yet. In particular, the characteristic plasma current (I p ) spike observed during disruptions is typically much smaller in simulations than in experiments [18][19] [20].

This paper presents JOREK simulations of an argon-MGI-triggered disruption in JET and their comparison to experimental data (for a recent overview on disruptions and MGI experiments in JET, see [21]). This work was initiated at a time when it was not yet clear that SPI would be selected for ITER. Future work should clearly focus on SPI rather than MGI simulations. However, MGI and SPI probably involve common physics to a certain extent. In particular, when the gas or pellet contains impurities (which is the case in this paper), MHD activity is thought to be destabilized in both cases by the penetration of a cold front and by helical cooling inside magnetic islands. On the other hand, when the gas or pellet contains pure deuterium, MGI and SPI may behave differently due a much milder plasma cooling from deuterium and to the much deeper penetration of SPI [10].

It is observed that, depending on the setting of the argon source in the simulations, an I p spike comparable to experimental measurements may or may not be obtained. What seems to determine the outcome is whether or not the thermal drive of the m = 2, n = 1 (or 2/1 in short; m and n are the poloidal and toroidal mode numbers) tearing mode persists in time, even as magnetic stochasticity develops. The thermal drive of tearing modes, i.e. the fact that radiation cooling inside a magnetic island may drive the island to grow further, has been pointed out by Rebut and Hugon several decades ago, already in the context of disruptions [START_REF] Rebut | Plasma Physics and Controlled Nuclear Fusion Research[END_REF], and has recently been invoked to explain density limit disruptions [START_REF] Gates | [END_REF]. The mechanism is relatively easy to understand: if radiation is strong enough inside an island, the temperature will decay locally, causing an increase in resistivity and thus a decrease in current density inside the island. The latter, similarly to the decrease in bootstrap current for neoclassical tearing modes [24], will drive further island growth. The picture becomes however more complicated when magnetic stochasticity comes into play, as is the case towards the end of the pre-Thermal Quench (TQ) phase of MMI-triggered disruptions, according to 3D MHD simulations. Indeed, heat conduction and shear Alfvén wave propagation along stochastic field lines will tend to refill the island region with heat and current, such that it is not obvious a priori whether the thermal drive of tearing modes, and especially of the largely dominant 2/1 mode, may continue to operate. In this paper, after describing the experiment in Section 2 and the model in Section 3, we will present in Section 4 a set of JOREK simulations in which we vary the argon injection parameters, depositing more or less argon in the 2/1 island region. Depending on parameters, the thermal drive of the 2/1 mode will be seen to persist or not as stochasticity develops. Simulation results will be compared in some detail to experimental data in Section 5. The violent relaxation in one of the simulated cases, which appears to have a persistent thermal drive of the 2/1 mode and produces an I p spike comparable to the experimental one, will be analyzed in Section 6. Finally, Section 7 will summarize the results and discuss their implications as well as the remaining open questions.

Experiment

We simulate the argon-MGI-triggered disruption in JET pulse 85943. This is an Ohmic discharge with a toroidal field B t = 3 T, plasma current I p = 2 MA, central electron density n e0 = 2.1 × 10 19 m -3 and temperature T e0 = 3.3 keV from an experimental session on Runway Electrons (REs). Pure argon is injected into a healthy plasma from Disruption Mitigation Valve 1 (DMV1), which is located at the top of the machine, as shown in Fig. 1. DMV1 is pre-loaded at 33 bar, corresponding to an argon content of 21.4 bar.L, i.e. 5.3 × 10 23 argon atoms. It is important to note however that, according to estimates based on gas dynamics theory [30], only about 1% of this quantity has entered the vacuum vessel by the time the thermal quench happens.

Figure 2 shows the evolution of I p , P rad (the radiated power inferred from bolometry measurements, assuming axisymmetric radiation), n el (the line-integrated electron density) from 3 lines of sight of the interferometry system (the signal is not valid anymore after 6.3 ms), magnetic fluctuations dB/dt measured by a Mirnov coil, and the electron temperature T e measured by 4 channels of the fast ('KK3F') Electron Cyclotron Emission (ECE) radiometer. The presented ECE data was digitised with a 200 kHz rate and smoothed with a +/-100 µs triangular filter to reduce noise. The time origin here is the time at which the argon enters the vacuum vessel (22.5282 s in 'JET absolute time'), calculated from the DMV1 opening time (22.5239 s) and the theoretical gas propagation time from DMV1 to the vessel (4.3 ms) [30] [31]. The pre-TQ phase lasts about 6 ms and is characterized by a gradual increase of n el and decrease of I p . The fast ECE data strongly suggests that a cold front penetrates during this phase: T e is seen to rapidly collapse for each channel, at time which increases from edge to central channels. The TQ is indicated by the collapse of T e from the central ECE channel. This follows a sharp increase in P rad up to about 500 MW within 1 ms. A burst of magnetic fluctuations is visible on the Mirnov coil signal for a period of about 1 ms after the TQ, a period during which P rad keeps rising fast, eventually peaking at 2.3 GW. I p also peaks, approximately at the same time as P rad . During the ensuing Current Quench (CQ) phase, a Runaway Electron (RE) beam is generated, resulting in a short I p plateau at about 1 MA around 12 ms.

Model

In order to simulate this experiment, we use the JOREK 3D non-linear magnetohydrodynamics (MHD) code, which is extensively described in [10] [25] and references therein. The specific model used in this work is the same as described in [9]. It is a reduced MHD model which includes an equation for the impurity density (argon here). The latter is summed over all charge states. Impurities are transported by diffusion as well as by advection at the same velocity as the main ions. A Coronal Equilibrium (CE) assumption is used to get the impurity charge state distribution and the associated energy sinks and sources from radiation, ionization and recombination, using data from the open ADAS database [START_REF] Summers | Atomic data and analysis structure users manual[END_REF]. A Collisional-Radiative (CR) model evolving the charge state distribution according to ionization and recombination rates has been implemented recently in JOREK [START_REF] Hu | Collisional-radiative non-equilibrium impurity treatment for JOREK simulations[END_REF]. Comparisons between the CE and CR models show that the CE model typically predicts a somewhat slower thermal collapse during an MMI [START_REF] Hu | Collisional-radiative non-equilibrium impurity treatment for JOREK simulations[END_REF]. Future studies aiming at quantitative validation will use the CR model.

In terms of simulation parameters, for the resistivity we use η 2×η Sp (min[T e , 0.7]), where η Sp is the Spitzer resistivity and T e is the electron temperature in keV. This way, the resistivity is realistic for T e ≤ 0.7 keV, while very low resistivities, which are numerically challenging, are avoided. The factor 2 is a simplified way to account for the fact that only passing particles carry the current (instead of a constant factor, a profile would be more accurate; the value 2 used here corresponds roughly to the inverse of the passing fraction near the edge of the plasma, which is the most important region for the pre-TQ and early TQ phases). Ohmic heating is included. The perpendicular kinematic viscosity (which is independent of temperature here, in contrast with previous works on also used at moderate level in order to smooth structures finer than the grid resolution.

For the poloidal discretization, which is done via cubic Bézier finite elements [28][10], a flux-aligned grid extending into the Scrape-Off Layer (SOL) with moderate resolution (3242 elements) is used, and for the toroidal discretization, Fourier harmonics from n = 0 to n = 10 are included. The simulations are run in free-boundary mode, making use of the JOREK-STARWALL coupling [29][10], and including a rather realistic, although axisymmetric, model of the JET resistive wall. It can be noted, e.g. in Fig. 4 below, that the Last Closed Flux Surface (LCFS) in the simulations is defined by the boundary of the computation domain on the high field side midplane instead of the X-point. This is not expected to affect results dramatically but should be improved in future studies.

One of the main limitations of the model is that it does not treat gas dynamics and is therefore not able to simulate the propagation and penetration of neutral argon from DMV1 into the plasma (for modelling of this process, see [30][31] [32][33]). Instead, argon is deposited via a volumetric source term and transported by diffusion and convection at the same velocity as the main plasma. At the beginning of the simulations, the argon is always deposited in the SOL. More precisely, the argon source is a Gaussian centered on R = 3.5 m, Z = 1.1 m (indicated by the cross in Fig. 3 -note that in reality the major radius of the gas arrival is closer to 3.25 m) and on the toroidal position of DMV1, with a width of 6 cm poloidally and 1.4 rad toroidally. In some simulations, as will be detailed below, the argon source is moved into the 2/1 island region (precisely at R = 3.3 m, Z = 0.9 m) after a few ms. As in [17] [18], the argon injection rate follows the rate expected from gas dynamics theory [30] but we use a smaller DMV1 pressure P DM V than the experimental 33 bar as a simple way of accounting for the fact that only a fraction of the argon penetrates the plasma. In order to select P DM V for the simulations, a comparison to interferometry data via a synthetic dagnostic is used, as shown in Fig. 3. In these simulations, where the argon source is left in the SOL at all times, it appears that P DM V = 6 bar allows a rather good match to the central lines of sight (2 and 3), while P DM V = 3 bar is not sufficient. The rise in the measurement of line of sight 4, which runs across the SOL on the low field side, is clearly not reproduced in either simulation, which may be a consequence of the fact that the model does not account for the argon that enters the vessel but does not penetrate the plasma.

Impact of the amount of argon deposited in the 2/1 island region

In this section, we will compare the above P DM V = 6 bar simulation in which the argon source is left in the SOL at all times, which we call 'Case A', with simulations in which the source is suddenly moved into the 2/1 island region (precisely R = 3.3 m, Z = 0.9 m) at 4.47 ms, which we call 'Case B' and 'Case C'. In Case B, P DM V remains equal to 6 bar after the source movement, while in Case C, it is increased to 12 bar. We will also analyze 'Case D', in which P DM V = 3 bar at all times and which has a source movement into the 2/1 island region at 5.02 ms. While Cases A, B and C may be directly compared to each other since their setups are clearly connected (they actually constitute a scan in the amount of material deposited in the 2/1 island region after 4.47 ms), Case D is not meant to be directly compared to the other cases. It is however included since, although less consistent with experimental data (e.g. interferometry or I p ) in the pre-TQ phase than the other cases (as will be shown below), it is more stable numerically and allows investigating the current profile relaxation and the I p spike after the TQ, as will be done in Section 6. The characteristics of these four cases are summarized in Table 1. Fig. 4 shows the evolution of T e for Cases A, B and C in the early phase of the simulation, before the source is moved (the three cases are identical in this phase). Three times are shown: t = 0 (left), 2.78 (middle) and 4.33 ms (right). The penetration of a cold front is clearly visible. The overlayed Poincaré cross-sections indicate that a stochastic layer extending from the edge to roughly the q = 2 surface has formed by 4.33 ms.

The ensuing evolution of T e in Case A is shown in Fig. 5. A TQ takes place at about 6 ms, caused by the stochastization of the magnetic field over almost the whole plasma volume. This TQ is however incomplete, with T e remaining on the order of several hundred eV in the core. Also, in the late pre-TQ (t = 4.92 ms, left plot) and early TQ (t = 5.95 ms, middle plot) phases, the 2/1 island region is not cooled to very Poloidal cross-section of the electron temperature, with a Poincaré crosssection (white dots)and the q = 2 and LCFS (calculated for the axisymmetric part of the magnetic field, green lines) overlayed, in the injection plane at t = 0 (left), 2.78 (middle) and 4.33 ms (right) for Cases A, B and C, which are identical in this phase (see Table 1). The color scale is saturated at 0.5 keV so as to highlight the dynamics in the cold regions.

Figure 5. Poloidal cross-section of the electron temperature, with a Poincaré crosssection (white dots) and the q = 2 and LCFS (calculated for the axisymmetric part of the magnetic field, green lines) overlayed, in the injection plane at t = 4.92 (left), 5.95 (middle) and 6.39 ms (right) for Case A (see Table 1). The color scale is saturated at 0.5 keV so as to highlight the dynamics in the cold regions. low temperatures.

Fig. 6 shows similar plots for Case B, in which the source is moved into the 2/1 island region at 4.47 ms. It is important to stress that moving the source in such a way has some justification. Indeed, as mentioned above, our model does not include gas dynamics. Simulations including gas dynamics [33] [31] find that when the gas front first reaches the plasma edge, only a small fraction of the gas penetrates the plasma. Poloidal cross-section of the electron temperature, with a Poincaré crosssection (white dots) and the q = 2 and LCFS (calculated for the axisymmetric part of the magnetic field, green lines) overlayed, in the injection plane at t = 4.62 (left), 4.91 (middle) and 5.09 ms (right) for Case B (see Table 1). The color scale is saturated at 0.5 keV so as to highlight the dynamics in the cold regions. This is the reason why we position the argon source in the SOL at the beginning of the simulation. The level of argon penetration in this phase is determined by P DM V and the argon diffusion coefficient, which we set so as to match interferometry measurements. However, as the cold front penetrates into the plasma, it is likely that the cloud of neutral argon outside the plasma expands radially. Hence, moving the argon source inward as the cold front penetrates appears reasonable, justifying to some extent our displacement of the source into the 2/1 island region at 4.47 ms. Comparing Figs. 5 and6, it can be seen that when the argon source is moved into the 2/1 island region, T e remains much lower there as global stochasticity develops. This is even more pronounced in Case C, where P DM V is increased to 12 bar after the source is moved, as shown in Fig. 7.

Fig. 8 compares the evolution of the magnetic energy in the n = 1 component (which is dominated by the 2/1 mode) for Cases A to D. It can be seen that the n = 1 magnetic energy shoots up in Cases B and C compared to Case A almost as soon as the source is moved. At about 5 ms, however, Case B 'rolls over'. On the other hand, the mode continues growing fast for Case C. Case D (which we recall has P DM V = 3 bar and a later source movement than in the Cases B and C, at 5.02 instead of 4.47 ms) shows a similar evolution to Case C but delayed in time.

Differences in the strength of the non-axisymmetric thermal drive are likely to at least partly explain these results. However, differences in the evolution of the current density profile may also have an effect. In order to help distinguish the two mechanisms, we ran simulations in which only the axisymmetric part of the resistivity η is taken into account, removing the non-axisymmetric thermal drive mechanism. Two such simulations with P DM V = 6 bar are shown in Fig. 9, one in which the source is left in the SOL at all times, equivalent to Case A (plain line), and one in which it is moved Poloidal cross-section of the electron temperature, with a Poincaré crosssection (white dots) and the q = 2 and LCFS (calculated for the axisymmetric part of the magnetic field, green lines) overlayed, in the injection plane at t = 4.62 (left), 4.80 (middle) and 5.09 ms (right) for Case C (see Table 1). The color scale is saturated at 0.5 keV so as to highlight the dynamics in the cold regions. 1).

into the 2/1 region at 5.95 ms, equivalent to Case B, although the source is moved later here (dashed line). The behaviour is similar to Cases A and B: the n = 1 energy shoots up in the second simulation after the source is moved and rolls over a few hundred microseconds later. This suggests a likely role of the current profile evolution in the 'shoot up' of the n = 1 energy for Cases B and C right after the source movement at 4.47 ms. Coming back to Fig. 8, one can however notice a striking bifurcation between Cases B and C shortly before 5 ms. With the support of Figs. 6 and7, we interpret this bifurcation as due to the persistence or not of the thermal collapse in the 2/1 island, and thus of the thermal drive. We note that Case D behaves similarly to Case C with about 1 ms of delay. A persistent thermal drive is thus likely to be at play also in Case D, although further work is needed to confirm this. Case D will be studied in more detail in Section 6.

Detailed comparison to experimental data

Fig. 10 compares the evolution of I p in the experiment and in Cases A to D. A first observation which can be made is that the pre-TQ evolution of I p in Cases A to C is quite similar to the experiment. Case D has a slower pre-TQ I p decay due to the lower P DM V . A second observation is that all simulations produce an I p spike, but with a clearly different magnitude depending on the case. Cases A and B have a much smaller spike than the experiment. Unfortunately, Case C runs into a numerical issue which prevents from observing the entire I p spike. Most interestingly, Case D produces an I p spike very similar to the experiment, both in magnitude and timescale. It is not clear why Case D has a much larger I p spike than the other cases. In fact, one may have expected the opposite since Case D has a lower P DM V . It should however be recalled that, while Cases A to C constitute a 'clean' scan in the amount of argon deposited in the 2/1 island region after 4.47 ms, Case D cannot be viewed as part of this scan. More specifically, the state of the plasma at the moment when the argon source is moved into the 2/1 island is not the same in Case D as in the other cases (e.g. Case D has slightly broader temperature and current profiles), and this must somehow explain the observed behaviour. In any case, we note that the larger n = 1 magnetic energy reached by 1).

Case D compared to Cases A and B (see Fig. 8) is consistent with its larger I p spike. We also note that Case C reaches an even larger n = 1 magnetic energy than Case D. Thus, one may think that Case C would have produced a similar or larger I p spike than Case D, had it not run into a numerical issue. When comparing only Cases A to C, a consistent trend of a larger I p spike (and n = 1 magnetic energy, see Fig. 8) with an increasing amount of argon deposited in the 2/1 island region after 4.47 ms seems to exist, although this is again slightly speculative since Case C runs into a numerical issue before completing its I p spike.

As will be discussed below, one non-realistic aspect of the simulations is that they do not capture the global radiative collapse at the end of the TQ which leads to the fast CQ. Thus, there is very little resistive I p decay during the I p spike in the simulations, while in the experiment, one may estimate (assuming that the time of the thermal collapse corresponds to the peak radiated power and that the resistive decay rate can be extrapolated from the CQ phase) that about 30 kA are resistively lost during the I p spike phase. This represents about 20% of the I p spike height, which is not negligible but not sufficient to strongly affect the comparison between simulations and experiments made in Fig. 10.

A final note related to the evolution of I p is that supra-thermal electrons during the TQ have been observed recently on DIII-D in sufficient amounts to carry potentially a significant fraction of I p [34]. These supra-thermals are not included in the present model. Their possible influence on the I p spike or more generally on the TQ and current profile relaxation is presently unclear and will have to be studied.

Another important point of comparison with the experiment is the radiated power, P rad . Fig. 11 compares P rad in Cases A to D with the experiment. A sharp rise is visible 1).

in Cases B and C when the argon source is moved into the 2/1 island region. As could already be guessed from Fig. 10, it is evident that the source is moved too early in Cases B and C. Moving the source later, like in Case D, provides a much better match in the evolution of both I p and P rad , with the rise of both quantities happening only a few hundred microseconds earlier in Case D compared to the experiment. What is clearly not reproduced by any of the is the continuing rise in P rad after 6 ms. This is probably due to a lack of global radiative collapse in the simulations, stemming from an insufficient argon influx. Unfortunately, as illustrated by Case C, it is numerically challenging to increase the argon source in the simulations.

In addition to the global P rad , we also compare the spatial structure of the radiation pattern. Fig. 12 shows a schematic view of JET from the top, indicating the position of DMV1 and of some relevant diagnostics, including the two 'KB5' bolometers (horizontal, KB5H, and vertical, KB5V). The latter are located at two distinct toroidal positions, both quite distant from DMV1. Their lines of sight are shown in Fig. 14.

Fig. 13 compares experimental (top) and synthetic (bottom) measurements for KB5H (left) and KB5V (right). The simulation used here is Case A, but the pattern is rather similar in all simulations. A striking geometrical resemblance can be seen in the experimental and synthetic patterns. In both cases, two maxima are visible for KB5H around 160 and 200 degrees in the pre-TQ phase, which gradually get closer to each other. For KB5V, a single maximum, near 270 degrees and moving towards smaller angles, is visible.

Fig. 14 shows cross-sections of the radiated power density in the simulation at 6.24 ms in the poloidal planes of KB5H (left) and KB5V (right), with corresponding lines of sight overlayed. These cross-sections help interpret the patterns of Fig. 13: they result from radiation being large at the edge of the 2/1 island region, where T e maximizes the argon radiative cooling rate. Note that, even though the parallel flow is somewhat damped and smoothed by parallel hyper-viscosity in these simulations (for numerical stability reasons), a non-negligible parallel flow still exists which makes the argon plume expand, explaining the radiation away from the injection point.

We then compare the experimental and synthetic n = 1 locked mode signal in Fig. 15, both in terms of amplitude (left) and phase (right). The synthetic signal is obtained via saddle coils, exactly like in the experiment (synthetic saddle coils with the exact geometry are included in the simulation). It can be seen that the amplitude of the locked mode signal has the right order of magnitude for all simulations. The amplitude is not extremely different between simulations, e.g. between Cases A and D, in spite of their large difference in terms of n = 1 mode energy (see Fig. 8) and I p spike (see Fig. 10). This may be related to a low-pass filter effect due to the finite wall penetration time (the saddle coils are just behind the wall). The phase of the locked mode signal is not well reproduced in the pre-TQ phase, but quite well at the TQ, i.e. at about 6 ms: 4 rad in simulations vs. 3 rad in the experiment. This indicates that, as it grows to a very large amplitude, the mode locks at approximately the same position in the experiment and simulations. The O-point of the 2/1 island coincides with the position where the gas is deposited, which is consistent with a thermal drive (i.e. a drive from helical cooling) and is also well known experimentally in cases with low toroidal rotation and weak error field [42][43]. This is also consistent with the good geometrical match observed in the radiation pattern (Fig. 13).

Finally, in Fig. 16, we compare experimental and synthetic T e profiles as measured by the fast ECE system. The experimental data is believed not to be affected by cut-off issues, based on estimates using n e profiles from the JOREK simulations, which we recall match interferometry measurements (Fig. 3). It is also not believed to be affected by non-thermal electrons, which usually produce a fast rise of the signal measured by some ECE channels. The synthetic data in Fig. 16 are from Case A. Several observations can be made on this figure. First of all, the central plasma is much hotter in the simulation than in the experiment. This is due to a spurious heating at the beginning of the simulation, which comes from the large particle diffusion used here. Indeed, since particle diffusion is implemented in such a way that heat does not diffuse with the particles, the diffusion of the initially peaked density profile causes the temperature profile to peak. A second observation is that the cold front penetrates faster in the simulation than in the experiment, and that the flattening of the T e profile around the q = 2 surface (near R = 3.4 m) happens earlier and in a more gradual way in the simulation than in the experiment: in the simulation profiles, one can see the flattened region growing between 4.80, 5.30 and 5.80 ms, while in the experiment, profiles at 4.80 and 5.30 ms are almost identical, with a substantial flattening appearing at 5.80 ms, directly followed by a complete T e collapse at 6.30 ms. It thus seems that the 2/1 mode appears later but grows faster in the experiment than in the simulation. Altogether, these instructive comparisons with experimental measurements suggest that simulations capture well the fact that a large 2/1 mode is destabilized by the MGI, with the 2/1 island O-point corresponding to the gas deposition region. However, the detailed pre-TQ dynamics are not well captured. In particular, the phase of the mode and the T e profiles from ECE are not reproduced. This may be the object of future studies. 6. Relaxation resulting from a persistent thermal drive of the 2/1 tearing mode

We will now look in detail at Case D, which produces an I p spike comparable to the experimental one, as seen in Fig. 10. Fig. 17 shows a series of poloidal cross-sections of the toroidal current density j φ (left) and its axisymmetric part j φ,n=0 (middle), as well as the electron temperature T e (right) at times 5.26, 5.72, 6.05 and 6.54 ms (top to bottom). Fig. 18 shows, for the same times, profiles as a function of the normalized n = 0 poloidal magnetic flux, ψ N ≡ ψ n=0 -ψ n=0,a ψ n=0,b -ψ n=0,a (where ψ n=0,a and ψ n=0,b are the n = 0 flux on axis and at the LCFS), of the flux surface averaged toroidal current density j φ (top left), safety factor q (top right), electron temperature T e (bottom left) and mass density (including all species) ρ (bottom right). Fig. 19 shows, still at the same times, profiles at the midplane, as a function of the major radius, of the n = 0 component of the toroidal current density j φ,n=0 (left) and electron temperature T e,n=0 (right). A substantial growth of the 2/1 mode between 5.26 and 5.72 ms, accompanied by the stochastization of the magnetic field over most of the plasma volume, is visible by comparing the first two rows of Fig. 17. This corresponds to an increase in the n = 1 magnetic energy by about one order of magnitude (see Fig. 8). Magnetic stochastization leads to a substantial flattening of T e between these two times (see Figs. 18 and19). However, the current density, safety factor, and mass density profiles have virtually not changed by 5.72 ms. The fact that T e responds faster to magnetic stochasticity than j φ or ρ is indeed expected and comes from the fast thermal motion of electrons along field lines. On the other hand, between 5.72 and 6.05 ms, a further growth of the 2/1 mode leads to a flattening of the j φ , q and ρ profiles in a large part of the plasma (ψ N ≤ 0.6, roughly). One can see in particular that q flattens to a value slightly above 2 in that region. As it does so, the q = 2 surface (indicated by the inner green line in Fig. 17) moves towards the plasma center. The j φ cross-section shows a thin current ribbon at the plasma center (and consistently the current profiles of Figs. 18 and 19 show a central peak), which we interpret as a consequence of the 2/1 'ghost island' 'running into itself' in the center and of the associated reconnection. These observations are reminiscent of the results of Drake and Kleva [36]. However, in the latter work, the 2/1 island runs into itself by growing so much that it encompasses virtually the entire plasma cross-section, which requires a central safety factor q 0 ≥ 1.5. In contrast, in the present simulation, q 0 is close to 1 and it appears that the 2/1 island runs into itself by moving towards the center along with the q = 2 surface as a result of the current profile evolution. A precise understanding of these dynamics remains to be established.

One can see in Figs. 17, 18 and 19 that the j φ profile has substantially flattened at 6.05 ms in the relaxed region. As shown in Fig. 20, this process takes place at almost constant global magnetic helicity, which is expected from theory for fast MHD relaxations [37][38][39] (we note that T e 0.2 keV during this phase, implying little resistive dissipation of helicity). This leads to an increase of the total current in the relaxed region. However, this increase is compensated by the induction of a negative skin current at the edge (clearly visible on the j φ,n=0 cross-section in Fig. 17), such that the total I p in the simulation domain does not change immediately (see Fig. 10). Consistently with the picture described in [START_REF] Biskamp | Nonlinear Magnetohydrodynamics[END_REF], the I p spike seems related to the decay of the negative skin current, which takes place between 6.05 and 6.54 ms (compared the last two rows in Fig. 17).

It can be noticed on Fig. 18 (bottom right) that the relaxation process also leads to a flattening of the mass density in the plasma core. This is associated to a strong mixing by an essentially macroscopic E × B flow associated with the 2/1 mode.

It is interesting to look into field line transport properties in this simulation. These may have important consequences regarding RE generation and may also help understand the current profile relaxation process [40][38] [39]. A simple analysis is done here: we initialize 200 field lines at R = 3.2 m, Z = 0.2 m (which corresponds to ψ N 0.22), evenly spread toroidally, and we track them for 200 toroidal turns. We then represent in Fig. 21, as a function of the number of toroidal turns, the minimal and maximal ψ N over all field lines (left), as well as the number of field lines remaining within ψ N ≤ 1 (right) at the same times as in Figs. 17, 18 and 19. It can be seen that as the 2/1 mode grows, field lines spread over a large part of the plasma (say ψ N ≤ 0.8) within a handful of turns. However, it takes roughly an order of magnitude larger number of turns for field lines to be lost. We attribute this difference to the fact that the 2/1 mode amplitude peaks in the core (see Fig. 9 in [40], which is from a similar JOREK simulation) and to the larger magnetic shear near the edge. The median number of turns for field lines to be lost decreases strongly with time: it is about 200 at 5.72 ms, 100 at 6.05 ms, and 20 at 6.54 ms. The reasons for this evolution are not entirely clear but may be related to the evolution of the modes amplitude and of the q profile (see Fig. 18).

Summary and outlook

The results presented above suggest that the argon MGI-triggered disruption in JET pulse 85943 is dominated by a 2/1 mode which reaches a very large amplitude as a consequence of the thermal drive persisting even as magnetic stochasticity develops. This triggers a relaxation of the current density profile in (roughly) the q ≤ 2 region, whereby q is brought to a value slightly above 2 in this region. I p spikes shortly after this relaxation, as the negative skin current induced at the edge by the relaxation decays. Although simulations do not match experimental measurements in detail, especially during the pre-TQ phase, a qualitative and sometimes quantitative match is obtained regarding several aspects: pre-TQ line integrated density evolution, late pre-TQ P rad and its spatial distribution, locked mode amplitude and phase at the TQ, pre-TQ I p evolution, I p spike. None of the simulations however produces a good agreement for all of these together, but this might be improved by moving the argon source later in simulations starting with P DM V = 6 bar. Also, none of the simulations presented here matches the experimental rise in P rad after the TQ. As mentioned above, this probably requires increasing the argon source at some point, which will be tried in the future (first attempts however ran into numerical issues).

We speculate that the dynamics observed in these simulations take place to some extent in a large number of disruptions triggered by MMI (whether MGI or SPI). However, at least three points remain to be clarified: 1) the role of toroidal rotation, 2) the role of pre-existing MHD modes and 3) the role of the 1/1 mode. Let us briefly discuss these. In the present paper, the target plasma (which is Ohmic) has a relatively slow toroidal rotation and does not suffer from pre-existing MHD modes. These conditions are observed experimentally to favor a direct influence of the MGI on the phase of the 2/1 mode, whereby the O-point of the 2/1 island tends to lock to the gas deposition region, as discussed in [42][43]. This is consistent with the strong thermal drive of the 2/1 mode found in the above simulations. However, with faster toroidal rotation (e.g. in the presence of substantial neutral beam heating), the 2/1 mode is often observed to rotate during the whole pre-TQ phase, possibly locking only at the beginning of the TQ [START_REF] Lehnen | [END_REF] [43]. On the other hand, a pre-existing locked mode may be the dominating factor determining the phase of the 2/1 mode, limiting the influence of the injection location [START_REF] Lehnen | [END_REF]. The 1/1 mode is not present in the above simulations since q 0 is just above 1, but former simulations with q 0 slightly below 1 gave a qualitatively similar behaviour (see slides 19 and 20 in [41]). This would suggest a small influence of the 1/1 mode. However, this mode has been found to play an important role in other simulations (see e.g. [9][16]) as well as in experiments (in relation to cold bubble observations [START_REF] Lehnen | [END_REF]). We speculate that the influence of the 1/1 mode depends on the radius of the q = 1 surface.

Finally, simulations take roughly one month on a High Performance Computing (HPC) cluster, which makes progress rather slow. Furthermore, it is numerically challenging to run through the TQ, which constitutes a bottleneck and will be the object of future efforts. In particular, we plan to assess how much numerical dissipation coefficients can be increased without making simulations unrealistic.

Figure 1 .

 1 Figure 1. Cross-section of the JET vacuum vessel showing the geometry of DMV1 (courtesy of S. Jachmich).

Figure 2 .

 2 Figure 2. Overview of the disruption phase of JET pulse 85943. From top to bottom: plasma current I p , radiated power P rad , line-integrated electron density n el from 3 lines of sight of the interferometry system (the signal is not valid anymore after 6.3 ms), magnetic fluctuations dB/dt (given here in terms of Mirnov coil voltage), and electron temperature T e for 4 channels of the fast ECE radiometer (the radii correspond to the mapped equilibrium at 22.3852 s in 'JET absolute time').

Figure 3 .

 3 Figure 3. Left: Comparison between experimental (plain lines) and synthetic interferometry data from two JOREK simulations with DMV1 pressures of 6 (dashed) and 3 bar (dash-dotted). Right: Cross-section of the electron density from the 6 bar simulation showing the interferometry lines with the same color code as in the left figure.

Figure 4 .

 4 Figure 4.Poloidal cross-section of the electron temperature, with a Poincaré crosssection (white dots)and the q = 2 and LCFS (calculated for the axisymmetric part of the magnetic field, green lines) overlayed, in the injection plane at t = 0 (left), 2.78 (middle) and 4.33 ms (right) for Cases A, B and C, which are identical in this phase (see Table1). The color scale is saturated at 0.5 keV so as to highlight the dynamics in the cold regions.

Figure 6 .

 6 Figure 6. Poloidal cross-section of the electron temperature, with a Poincaré crosssection (white dots) and the q = 2 and LCFS (calculated for the axisymmetric part of the magnetic field, green lines) overlayed, in the injection plane at t = 4.62 (left), 4.91 (middle) and 5.09 ms (right) for Case B (see Table1). The color scale is saturated at 0.5 keV so as to highlight the dynamics in the cold regions.

Figure 7 .

 7 Figure7. Poloidal cross-section of the electron temperature, with a Poincaré crosssection (white dots) and the q = 2 and LCFS (calculated for the axisymmetric part of the magnetic field, green lines) overlayed, in the injection plane at t = 4.62 (left), 4.80 (middle) and 5.09 ms (right) for Case C (see Table1). The color scale is saturated at 0.5 keV so as to highlight the dynamics in the cold regions.

Figure 8 .

 8 Figure 8. Time evolution of the magnetic energy in the n = 1 component for a set of simulations with different argon source settings (see text and Table1).

Figure 9 .

 9 Figure 9. Time evolution of the magnetic energy in the n = 1 component for two simulations (equivalent to Cases A and B from Fig. 8) in which only the axisymmetric part of the resistivity η is used.

Figure 10 .

 10 Figure 10. Evolution of the plasma current in the experiment and in a set of simulations with different argon source settings (see text and Table1).

Figure 11 .

 11 Figure 11. Evolution of the radiated power inferred from bolometry measurements in the experiment and in a set of simulations with different argon source settings (see text and Table1).

Figure 12 .

 12 Figure 12. Schematic view of JET from the top, indicating octant numbers and the position of DMV1 and of interferometry, bolometry and ECE systems.

Figure 13 .

 13 Figure 13. Experimental (top row) and synthetic (from Case A, bottom row) bolometry measurements for the horizontal (KB5H, left column) and vertical (KB5V, right column) bolometers. Please note the different color scales.

Figure 14 .

 14 Figure 14. Poloidal cross-section of the radiated power per unit volume from Case A at 6.24 ms in the plane of the horizontal (KB5H, left) and vertical (KB5V, right) bolometers, with lines of sight shown in white and extreme angles indicated for easier comparison with Fig. 14.

Figure 15 .

 15 Figure 15. Experimental and simulated (for a set of simulations with different argon source settings -see text andTable 1 locked mode signal amplitude (left) and toroidal phase (right).

Figure 16 .

 16 Figure 16. Comparison between experiment (thin lines with crosses) and synthetic (thick lines) electron temperature from fast ECE. The synthetic data is from Case A.

Figure 17 .Figure 18 .

 1718 Figure 17. Poloidal cross-sections in the MGI plane of the toroidal current density j φ (left), the axisymmetric part of the toroidal current density j φ,n=0 (middle), and the electron temperature T e (right), with a Poincaré cross-section (white dots) and the q = 2 and LCFS (calculated for the axisymmetric part of the magnetic field, green lines)

Figure 19 .

 19 Figure 19. Profiles vs. the major radius of the n = 0 component at the midplane of the toroidal current density (left) and the electron temperature (right) for Case D at 5.26, 5.72, 6.05 and 6.54 ms (same times as in Figs. 17 and 18).

Figure 20 .

 20 Figure 20. Evolution of the experimental and simulated plasma current as well as simulated plasma internal inductance and global magnetic helicity for Case D.

Figure 21 .

 21 Figure 21. Field line statistics for Case D, based on a population of 200 field lines initialized at R = 3.2 m, Z = 0.2 m (which corresponds to ψ N 0.22) at evenly spread toroidal positions.

Table 1 .

 1 Characteristics of the simulation cases.

Table 1

 1

locked mode signal amplitude (left) and toroidal phase (right).

MMI-triggered disruptions with JOREK[17][18][6][7][8][9]) and heat diffusivity are of order 3 m 2 /s and 2 m 2 /s respectively, which might be considered as representative of underlying turbulent transport, although the characteristics of the latter during an MGI are unknown. A Spitzer-Härm parallel heat conductivity is used. The most unrealistic aspects of the model are probably the particle diffusivity, set to the large value of 30 m 2 /s for both main ions and argon, and the fact that we artificially damp the parallel flow via a large hyper-parallel-viscosity. These choices have been made in order to avoid numerical instabilities. Hyper-resistivity and hyper-perpendicular-viscosity are

8. Acknowledgements Parts of this work have been carried out within the framework of the EUROfusion Consortium and have received funding from the Euratom Research and Training Programme 2014-2018 and 2019-2020 under Grant Agreement No. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. Part of this work was performed using HPC resources from GENCI-TGCC (Grant A8-gen11358 ). Part of this work was carried out using the Marconi-Fusion supercomputer. Part of this work is supported by the National Natural Science Foundation of China under Grant No. 11905004.

We thank Michael Lehnen, Stefan Jachmich and Uron Kruezi for providing information on the DMV, Dylan Brennan and Patrick Maget for providing some of the references, and Cédric Reux for useful discussions.