Development of full electromagnetic plasma burn-through model and validation in MAST

This paper describes the improvement of the electromagnetic plasma burn-through model. Full circuit equations describing the currents in solenoid, poloidal field coils, and toroidally conducting passive structures have been integrated into the differential equation system of the plasma energy and particle balances in DYON. This enables consistent calculation of the time-evolving loop voltage at a plasma position only using operation signals in a control room, which are current (or voltage) waveforms in solenoid and poloidal field coils and prefill gas pressure. The synthetic flux loop data calculated in the modelling agrees well with the measurement in MAST, confirming the validity of the loop voltage calculation. The electromagnetic modelling also enables calculation of 2D time-evolving poloidal magnetic flux map, thereby modelling the plasma volume evolution during the plasma break-down and burn-through phase. Only using the control room operation signals used in 34 ohmic start-up discharges with the direct induction start-up scenario in MAST, the electromagnetic plasma burn-through modelling has reproduced the time-evolution of plasma current, electron density and temperature, and plasma volume, showing a reasonable level of agreement with experimental measurement.


Introduction
Tokamak start-up consists of three phases: plasma breakdown, plasma burn-through, and plasma current ramp-up. The induced high loop voltage or injected electromagnetic wave enables plasma break-down, generating a plasma at a * Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. low degree of ionization in a prefilled hydrogen isotope gas. Since the degree of ionization is low in this phase the electron-neutral collisions are dominant [1]. This justifies the Townsend break-down theory, which calculates the minimum toroidal electric field required for plasma break-down E Townsend at a given prefill gas pressure p and an effective length of open magnetic field lines L f i.e. 1 = α(p, E Townsend )L f [2] where α is the mean number of collisional ionizations by an electron in 1 m (the explicit equation indicated in (13)). As the degree of ionization increases, the energy loss rate of electrons and ions rapidly rises due to the increase in atomic reactions such as ionization, radiation, and charge exchange. This impedes further increase in T e , and high ohmic heating power or electron cyclotron heating (ECH) is required to overcome the energy loss until achieving sufficient ionization. Full ionization of deuterium could be achieved at T e less than 5 (eV) [3], but fully stripping low Z impurities e.g. C 6+ , or O 8+ (i.e. impurity burn-through) requires a few tens of eV [4]. After passing the peak energy loss in the plasma burn-through phase, the reduced electron energy loss enables the increase in T e , thereby increasing plasma current with decreasing plasma resistance. Since the loop voltage or ECH power needed for sufficient plasma burn-through is much higher than that for plasma break-down, the main hurdle for successful tokamak start-up is in the plasma burn-through phase [5].
Operation scenarios for tokamak start-up are usually sought with the trial-and-error method when a new device is built. Once a few operation recipes are found in the device, the scenario is usually not further optimized. So far, this approach has not had a serious problem in present devices, although it often made a few months of delay to starting scientific experiments. However, it is likely that the approach is not enough in future devices. It is foreseen there will be a challenge in plasma initiation in ITER, which is a large superconducting device. The maximum toroidal electric field available in ITER is about 0.3 (V m −1 ) [5], which is much lower than the value typically used in present devices e.g. 1 V m −1 in JET. In addition, the large L/R time τ L/R 5 resulting from the large thick vacuum vessel in ITER would make the plasma control complicated [6]. Considering about a decade of delays of ITER construction, immediate successful generation of ITER first plasma without any operational issue is of crucial importance. For this reason, it is now vital to optimize the ITER start-up scenario with reliable plasma burn-through modelling.
Plasma burn-through modelling is also important to design a new device such as spherical tokamak for energy production (STEP) [7]. In a spherical tokamak, due to the limited space and the high neutron irradiation in the centre column it is beneficial to have a thin or even no solenoid. However, this would significantly increase the uncertainty of plasma burnthrough. In order to design a new spherical tokamak, it is essential to assess with numerical modelling the feasibility of plasma initiation and current ramp-up with a limited or zero loop voltage induced by the solenoid.
DYON is a plasma burn-through simulation code, developed and validated in JET with carbon wall [8] and ITER-like wall [9], and used for ITER prediction [10]. Consistency of DYON with other mainstream plasma burn-through simulation codes, SCENPLINT [11][12][13] and BKD0 [14][15][16], was verified by the code comparison conducted as a joint activity in International Tokamak Physics Activity (ITPA)-Integrated Operating Scenario (IOS) group. Brief history of development of the plasma burn-through codes is summarized in [17], which reports the code benchmark activity. One of the main findings 5 A characteristic time for eddy currents in the vessel structures to decay i.e.
in [17] is that the vessel eddy currents and the evolving plasma volume could significantly affect the simulation results.
In JET modelling with DYON [8,9], the eddy current was assumed only in the structure supporting the divertor (called MK2 structure) as it is the only toroidally conducting passive structure in JET. Two circuit equations (i.e. plasma current and the eddy current in MK2) were enough to simulate JET discharges reproducing measured experimental data in the modelling. However, such simplified circuit equations might not be sufficient in other devices including ITER which has more toroidally conducting structures. With this motivation, full circuit equations describing all toroidal conducting vessel structures, solenoid, and PF coils have been integrated in DYON. The upgraded DYON not only allows proper modelling of vessel eddy currents but also calculation of loop voltage at the plasma region with the input of control room operation scenarios such as power supply voltage or coil currents. Loop voltage calculation at the plasma region is of importance as it could differ from the measured value at the vacuum vessel due to the stray field, in particular, for tokamaks with no iron core e.g. ITER.
Plasma volume is an important parameter in the plasma burn-through models which solve the energy and particle balance in 0D i.e. calculation of volume-averaged T e and n e . Since there is no measured data available in the burn-through phase, in previous modelling the plasma volume was prescribed with assumed values, which allow reproducing experimental data such as D α emission. For example, the plasma volume in the code benchmark with JET data was given by interpolating the vacuum vessel volume at 0 s and the first EFIT plasma volume data [17]. This approach was useful to validate the plasma burn-through model, but would give a large uncertainty if the model is used to optimize operation scenario in future devices or to assess a new tokamak design. For predictive modelling, consistent calculation of plasma volume with a physics-based model is necessary. Hence, we have developed a plasma volume model in DYON. Implementation of the full circuit equations enables simulating the time-evolution of 2D poloidal magnetic flux map (i.e. ψ map) in the vacuum space. This allows assessing the Townsend break-down condition in each field line, and estimating the plasma volume. This paper will report the improvement and validation of an electromagnetic plasma burn-through model where the evolving plasma volume is consistently calculated. In section 2, the development of full circuit equations and the plasma volume model are introduced. Section 3 describes ohmic plasma initiation with the direct induction scenario in MAST, and compares simulation results against experimental data from a typical direct induction discharge. With the simulation settings used to reproduce the MAST discharge, the prediction capability of electromagnetic DYON simulation is assessed comparing the simulation results against 34 discharges, which were randomly selected among direct induction discharges in MAST. Section 4 provides the conclusion.

Full circuit equation
The RL circuit equation system describing electric currents in tokamaks consists of plasma current, active coil currents and induced currents in passive conducting structures, all of which are linked one another through mutual inductances. To calculate the loop voltage at the plasma region, the differential equation system of circuit equations should be solved.
We could first consider the simplest case which is a single coil loop without any other coils nearby. The consumption of the coil voltage applied by external power supply V c is decomposed into the resistive voltage (V res = R c I c ) and the inductive voltage (V ind = dψ c dt ), where I c and R c are the electric current and resistance in the coil loop, respectively. ψ c is the magnetic flux enclosed by the coil loop, and is calculated as ψ c = L c I c where L c is the self-inductance of the coil loop. If the coil geometry is fixed (i.e. dL c dt = 0), which is the case in experiment devices, the inductive voltage can be replaced with L c dI c dt : Now we consider multiple coil loops and toroidally conducting vessel structures (i.e. passive loops) nearby where magnetic fluxes are linked one another and the position of each loop is fixed. The system of circuit equations with n loops of coils and vessel structures have the following matrix form: where − → V cv is a n × 1 column vector of voltage applied by external power supply, to coil loops and toroidally conducting passive structure in the vacuum vessel (hence zero for the vector elements corresponding to the vessel components).
← → R cv is a n × n diagonal matrix of resistances in the coil loops and the vessel components.
− → I cv is a n × 1 column vector of currents in the coil loops and the vessel components.
← → M cv is a n × n symmetric matrix of inductances in coil loops and vessel components. In ← → M cv , non-diagonal elements and diagonal elements are mutual inductance and self-inductance, respectively.
In order to include the plasma current, (1) is extended by attaching the plasma current I p in the (n + 1)th row of the corresponding column vectors. In other words, − → V cv and − → I cv are extended to − − → V cv p and − → I cv p , which are (n + 1) × 1 column vectors. The plasma resistance R p is added in the (n + 1, n + 1) element of ← → R cv . The mutual inductances between I p and both coils and vessel components are attached in the (n + 1)th row and column of ← → M cv , and the plasma self-inductance is added in the (n + 1, n + 1) element. The extended mutual inductance matrix is ← − → M cv p , which is a (n + 1) × (n + 1) matrix.
← → R cv p , and ← − → M cv p can be explicitly written as (2) Compared to the circuit equation system without I p , one of the main differences in the extended circuit equation system is that the mutual inductance could evolve in time, as the position and shape of a plasma current could change in time.
This makes an additional term dt , the differential equation system describing all toroidal currents in a tokamak can be written: Note, R and Z are the radial and vertical position of each current element, respectively. All current loops are assumed to be toroidally symmetry. The elements of are all zero except the elements corresponding to the mutual inductances of I p with coils and vessel components i.e. the last row and column. This is because the positions of coils and vessel components are fixed (i.e. ∂R ∂t = 0 and ∂Z ∂t = 0), whereas I p position could change in time (i.e. ∂R ∂t = 0 and ∂Z ∂t = 0). Note, (n + 1, n + 1) element in is defined as zero. This is because the plasma self-inductance is not only a function of the radial position R and vertical position Z but also of the minor radius and elongation. Hence, the time derivative of plasma self-inductance is included as a separate dL p dt , and all the other elements are zero. L p is calculated with the evolving plasma parameters from the plasma volume model, which will be introduced in section 2.2.
If the power supply voltage applied to each coil is available to use as input data, the elements corresponding to coils in − − → V cv p (i.e. V c1 , V c2 , . . .) could be directly given while the other elements for vessel components and the plasma current are zero i.e. V v# = 0 and V plasma = 0, as there is no power supply connected to vessel components and I p . In this case, (3) can be solved for d − − → I cv p dt as it is. However, it is more common in tokamaks that coil current data is available to use, rather than power supply voltage data, as coil currents are easier to routinely measure. In the case of using coil currents as input data (i.e. − → I c and d − → I c dt ), the vessel currents and plasma current (i.e. d − → I vp dt ) can be calculated by solving a subset of the circuit equation matrix corresponding to the vessel components and plasma current. Note, the last row of (5) correspond to the circuit equation for I p , which can also be explicitly written: V loop is the loop voltage induced to I p by the coil currents and vessel eddy currents, and is calculated as − dψ cv dt , which is a n × 1 column vector of the magnetic flux induced by the time evolution of the currents in coils and vessel components T is the transpose of a n × 1 column vector − → M p , which corresponds to the mutual inductance of I p with coils and vessel components.
The plasma current position R and Z in the circuit equation could be calculated by solving radial and vertical force balance equations. However, such a force balance model is not yet available in the present version of DYON. In the present modelling R and Z are prescribed by input data. Note, when the field lines are closed, R corresponds to the magnetic axis position, and it could be different to the centre of the plasma cross-sectional area. The calculation of plasma self-inductance L p requires an assumption that the plasma shape is an ellipse, which is a function of the centre position of the plasma crosssectional area (i.e. R vol and r vol ), and elongation κ. R vol , r vol , and κ are calculated by the plasma volume model described in the next section.
Depending on the input data availability of coil voltage or current, either (3) or (5) is solved, together with the existing differential equation system in DYON. This enables consistent calculation of I p and currents in the coils and the vacuum components with time-evolving R p , which is determined by T e , Z eff , and the plasma volume model. Following is the summary of the differential equation system in DYON, which describes the plasma energy and particle balance. More details about DYON can be found in [8][9][10][17][18][19].
The radial and vertical magnetic fields and the induced loop voltage (and toroidal electric field) in each mesh grid cell are calculated by derivatives of ψ(R grid , Z grid ) with respect to R grid , Z grid , and time, respectively: B R (R grid , Z grid ) and B Z (R grid , Z grid ) in (10) and the toroidal magnetic field B φ (R grid , Z grid )(= , where R 0 is the major radius to the vacuum vessel centre) enables the field-linefollowing calculation. In other words, if R grid,k and Z grid,k are the radial and vertical position at kth point on an open magnetic field line, we can calculate the next mesh grid point after traveling a step length Δl along the field line i.e. R grid,k+1 and Z grid,k+1 Here, Δl is a constant length, and it should be sufficiently short considering the curvature of field lines. If not, the mesh grid points calculated with the field-line-following could deviate from the field line. However, too short Δl would require unnecessarily significant computation time 6 For open field lines, E is compared to the Townsend electric field E Townsend , which is the minimum electric field required for an electron to produce, at least, a single ionization event along an open field line, thereby triggering the Townsend break-down. E Townsend for a hydrogen isotope plasma is calculated as follows [5]  is considered to be infinite and the volume within the closed field lines are also counted. As will be shown later, the plasma volume is dominated by the open field lines at the beginning of breakdown, and as plasma current increases, it is gradually taken over by the volume of the closed field lines.
Note, the Townsend break-down assessment in each field line is to estimate the plasma volume, and the exponential increase of n e during the electron avalanche in the field lines are not modelled in DYON. Once E E Townsend in open field lines, DYON starts to use the calculated plasma volume to solve the energy and particle balance equation system in (8), and the initial plasma parameters are given by the predefined values in table 1 assuming the time for the electron avalanche is short enough to ignore. As will be shown by the prompt increase in n e in figure 6, the electron avalanche phase in the MAST discharges is completed within a few ms, and it is not necessary to calculate the time taken to reproduce experimental data in the plasma burn-through modelling. This should be due to the feature in a spherical torus where E φ /p (kV m −1 Torr −1 ) is high. It was reported that free acceleration of electrons in a spherical torus could shorten the time of electron avalanche phase [20]. However, in the case there is a significant delay in the electron avalanche for some reasons e.g. conventional tokamaks or space charge effects [21], the time for electron avalanche might also need to be considered in the model.
Once the plasma shape is determined by the plasma volume model above, the major radius of plasma volume R vol is calculated as the centre of mass, and the minor radius r vol is found by the distance from R vol to the first wall. The plasma elongation κ is estimated by κ = min A/(π * r 2 vol ), max(Z grid ) − min(Z grid ) /(2r vol ) , where A is the plasma cross-sectional area from the plasma volume model and Z grid is the vertical position of mesh grids occupied by the plasma. An ellipse approximated to the plasma shape is determined by R vol , r minor , and κ (see an example in figure 4). Using the ellipse, the plasma self-inductance for a conventional tokamak with a large aspect ratio [8] is calculated as: For a spherical torus with a small aspect ratio, which is the case in this paper, the plasma self-inductance formula is derived in [22,23] as: Consistent calculation of L open also improves the confinement time model with a better physics basis than the models used in present burn-through modelling codes (i.e. [17,24]), in which the evolving connection length during the transition from the open field to closed field lines was estimated by the increase in I p in the exponential term. In the new confinement time model, the parallel particle transport is calculated with the averaged length of open field lines L open , which is calculated with the field-line-following calculation in (11), and the ratio of the open field lines to the total field lines. The perpendicular particle transport is calculated with Bohm diffusion: #of open field lines #of total field lines where . (17)

Ohmic plasma initiation with direct induction scenario in MAST
There are two plasma initiation scenarios in MAST, namely merging compression and direct induction [20,25]. Figure 2 shows the solenoid and poloidal field coils in MAST. The former scenario was to trigger the plasma break-down near P3 coils by rapidly ramping up the currents in P3 coils and inducing a high loop voltage nearby. The P3 coil current decreases towards 0(A) after the peak value, and this makes the two plasma rings merged on the mid-plane. The merged plasma ring is further compressed with P4 and P5 currents to achieve plasma burn-through. After this, the magnetic flux swing in the solenoid induces the loop voltage to increase the plasma current. The latter scenario directly uses the magnetic flux swing in the solenoid and induces the loop voltage near the magnetic field null. This makes the plasma break-down and burn-through at the null, and also increases the plasma current. Although the merging compression scenario has a benefit of reducing the magnetic flux consumption in the solenoid (i.e. smaller Volt × seconds) and was well routinely used in most MAST discharges, the merging compression scenario would not be available to use in ITER or STEP as there will be no such in-vessel PF coils. For this reason, we have selected a MAST discharge operated with the direct induction scenario (#27512), and performed electromagnetic simulations using the upgraded DYON code. were given to − → I c and d − → I c dt in (5). The mutual inductance matrix and the resistance matrix corresponding to the MAST design were also given in (5), respectively. Solving the full circuit equations, the vessel eddy currents ( figure 3(b)), the plasma current ( figure 3( f )) and the evolution of 2D poloidal magnetic flux map (see figure 4) are consistently calculated. The evolving 2D poloidal magnetic flux map is used to calculate the evolving plasma volume with the model described in section 2.2. As can be seen in figure 4, due to the open magnetic field lines, the plasma volume is extended to the vessel wall at the beginning. The plasma volume shrinks to a smaller volume as the plasma current increases and the field lines are closed. The first plasma volume data in EFIT is available from 30 (ms). At around 30 (ms), the field lines are also fully closed in DYON modelling results (see figure 4). The calculated plasma volume agrees well with the EFIT data [26] (see figure 3(c)). Assuming toroidal symmetry of the plasma geometry, the modelled plasma volume is also used to find the plasma self-inductance and the plasma resistance (figures 3(d) and (e), respectively). These were used to calculate the plasma current evolution in the full circuit equations. The calculated plasma current agrees well with the measured value in figure 3( f ), implying the validity of the electromagnetic simulation.

Simulation of a MAST discharge
The simulation parameters used in the DYON modelling are summarised in table 1. As shown in figure 3(a), the rate of change in the solenoid and PF coil currents before t = −10 (ms) is small, and the induced vessel eddy currents at around −10 (ms) are small enough to ignore. DYON starts to solve the full circuit equation system from −10 (ms) with the initial condition that the vessel eddy currents in all passive structures are 0 (A). The positive loop voltage at the plasma region starts to increase from 0 (s). Since the D α peak appears at 3 (ms) in the photomultiplier tube measurement, the field-line-following calculation and the Townsend break-down model starts to calculate the plasma volume from 3 (ms), which are used to solve the differential equation system of the plasma energy and particle balance. In the MAST discharge modelling, the field-line-following calculation was performed for the magnetic field lines crossing 50 radial points on the mid-plane at every 1 (ms). The initial condition of the electron temperature T e (t = 0) was assumed as 3 (eV), and the initial degree of ionization γ(t = 0) was consistently calculated as 1% with γ(0) = 2 × 10 −3 T 1.5 e (eV) [8]. Both the prefilled D 2 molecular gas temperature T n and the initial deuterium ion temperature T i (t = 0) was assumed as the room temperature i.e. 0.026 (eV). This assumption was used to calculate the density of the prefilled deuterium atoms with the measured prefill gas pressure (Pa) at 0 s: It was found that use of somewhat lower measured gas pressure improves the agreement with experimental data in #27512. 75% of measured prefill gas was used as input. Since the gas pumping is not modelled, the injected external gas fuelling data, of which the peak value is about 1 × 10 22 (s −1 ), is too high to be used as it is. 10% gas fuelling efficiency was identified for DYON to reproduce the measured parameters in #27512. In addition, the measured gas pressure waveform in #27512 is about 10 ms behind the external gas fuelling waveform. Based on this, 10 ms delay has been taken into account in the input data of external gas fuelling (see figure 8).
The initial carbon content was assumed to be zero. Since the first wall material in MAST was made of carbon, carbon atom influx was calculated with the plasma-wall-interaction model validated against the JET data with the carbon wall in [8], where the chemical sputtering is defined as 0.03. To represent a tiny amount of the impurities remaining from previous discharges, 0.1% initial oxygen atom content in the prefill gas was assumed as an initial condition. The evolving density of all impurity charge stages was calculated in (8).
The loop voltage induction is a function of the rate of magnetic flux change in a toroidal loop, and thus depends on a loop position (i.e. radius and vertical position). Without an iron core, which is the case in most present devices including MAST, the loop voltage at a plasma position could be different from the value measured at the vacuum vessel, and it should be calculated by modelling. Figure 5(a) compares the loop voltage at the plasma position calculated by FIESTA [27] and DYON, using the same coil currents in figure 3(a). The So far, validation of plasma burn-through modelling has been only possible with a limited data set, as measured data are usually not available in such an early phase. It is because the main diagnostics system in present devices e.g. Thomson scattering (TS) is optimized for the flat-top phase where T e and n e are high. The Nd:YAG TS system [28] on MAST operates well for the start-up phase since it is designed to measure a low T e (from 5(eV)) and n e (from 1 × 10 18 (m −3 )) using its large collection solid angle. In MAST discharges, measured T e and n e are available from very early time e.g. 5(ms), and the profile data with high radial resolution covers the break-down and the burn-through phase i.e. 0-50(ms). This provides exceptional data to validate plasma burn-through models. Figures 6(a) and (b) shows the profiles of T e and n e measured by MAST TS system, where the line of sight is aligned with the mid-plane. It is worth noting that the first measured profiles of T e and n e at 5(ms) show that the plasma initiation takes place from the inboard side, which is qualitatively consistent with the modelled plasma volume in figure 4. Another useful point to note is that both T e and n e profiles are not much centrally peaked. One of the underlying assumptions in the mainstream plasma burn-through modeling codes (i.e. DYON, BKD0, and SCENPLINT) is that T e and n e are uniform across the radius in the burn-through phase. The measured TS profile data justifies this assumption. Figures 6(c) and (d) compares the modelled T e and n e in DYON against the time trace of the averaged TS data. The error bars in TS data have been calculated with the root-mean-square of the statistical errors in each radial point. The time evolution of T e and n e in DYON is within the experimental error bars, implying the discharge is well simulated in the modelling.
Plasma burn-through modelling could provide useful information to interpret the discharge in the burn-through phase. For example, figures 7(a) and (b) show how the electrons and ions lose their energy, respectively. At the beginning (i.e. ∼8 (ms)), electrons lose the energy for the deuterium radiation and the equilibration, which is the heating power to ions. At this time, ions mainly lose the energy for the charge exchange reactions. As can be seen in figure 7(c), deuterium is fully ionized at around 8(ms), and the deuterium radiation and the charge exchange energy loss are reduced. At 15-20(ms), as C 3+ becomes the dominant charge level, a much larger electron energy loss occurs due to the carbon radiation. The carbon radiation is significantly reduced as the dominant charge level further increases to C 4+ . After overcoming the carbon radiation, T e stably increases with the ohmic heating, which allows the plasma current ramp-up.

Statistical validation
To assess the prediction capability of the electromagnetic DYON modelling, 34 MAST discharges with direct induction ohmic plasma initiation were randomly selected. The same simulation parameters in table 1 were used in all the 34 DYON simulations, and the control room input data such as prefill gas pressure, external gas fuelling rate, and the coil current trajectories in the solenoid and PF coils were individually given from experimental data. Figures 8(a) and (b) show the measured gas pressure and the external gas fuelling rate in the 34 discharges. The prefill gas pressure at 0 s is in the range of 1.1-1.7(mPa), and has the peak value at around 15 (ms) due to the external gas fuelling, which was also used as input data. Figure 8(c) shows the loop voltage measured in the flux loop located at the inboard mid-plane (R = 0.18 (m) and Z = 0.015 (m)). The loop voltage waveform is almost identical until 20 (ms), indicating that the same direct induction scenario was used in all the 34 discharges. Considering that 34 discharges were randomly selected, this reveals that the plasma initiation scenario was only developed when the device was first explored after the construction, and has not been further optimized since one working recipe had been identified. This should be the case in most other present devices. For better validation of plasma burn-through models, dedicated plasma initiation experiments (e.g. operation parameter scanning for plasma initiation) are required. Among the 34 discharges, 11 discharges start feedback control at around 20 (ms) (the other discharges from 60 (ms)), and the scenarios used after 20 (ms) are varying. Since the modelling was performed until 50 (ms), the different loop voltage data in 20-50 (ms) are still useful for the statistical comparison. Figure 9 compares I p , n e , and T e calculated by electromagnetic DYON modelling against 34 MAST direct induction ohmic discharges. The colors of triangle symbols represent the time of comparison i.e. blue at 50 (ms), red at 40 (ms), cyan at 30 (ms), and green at 20 (ms). The plasma currents modelled in DYON agree very well with the measured values at all four comparison times (see figure 9(a)). The plasma currents at 20 (ms) are all closely gathered at around 120 (kA) in both calculated and measured. This is not surprising as the identical loop voltage data was used until 20 (ms) in all 34 discharges. The different increase in plasma currents are clearer with time. Although the measured plasma currents of 34 discharges at 50 (ms) are in a much wider range, they are still well reproduced in DYON, providing the confidence in I p ramp rate prediction in the burn-through phase.
The electron temperature calculated also reasonably agrees with TS measured electron temperature (see figure 9(c)), implying the plasma resistivity calculation should be in reasonable agreement as well. Compared to the very good agreement in the plasma currents, the scatter in figure 9(c) is larger. This could be due to the TS error bars as shown in figure 6(c) (e.g. 30 (eV) at 50 (ms)). Most electron densities in 34 discharges at the four comparison time are in the range of 5-10 × 10 18 (m −3 ). As can also be seen by the TS measurement in figure 6(d), the electron densities in 34 discharges rapidly increase until 20 (ms), but do not further increase after 20 (ms). A few discharges have measured densities higher than 12 × 10 18 (m −3 ). They correspond to the discharges with high external gas fuelling, which are also indicated by the increase in the gas pressure after 30 (ms) in figure 8(a) and (b). Since the external gas fuelling is one of the input data in the modelling, the increase in electron density for n e > 12 × 10 18 m −3 is qualitatively reproduced in the modelling as well, although they are about 30% undercalculated (see figure 9(b)). The lower n e prediction could be attributed to the constant effective gas fuelling coefficient, which is identically assumed as 10% in all discharges. Elaborated calculation of gas fuelling coefficients in individual discharges could improve the prediction capability.
Finally, figure 9(d) compares the calculated volume in DYON against the EFIT plasma volume. The first data of plasma volume in EFIT is available at 40 (ms), which is the time that the closed flux surfaces are fully formed. As shown in figure 4, the field lines at 40 (ms) are also fully closed in DYON modelling and the calculated volume agree very well with the EFIT data. For complete validation of the plasma volume model, measurement of plasma volume with open field lines is required. This would require a dedicated experiments or novel processing of fast camera images.

Conclusion
Following the lesson learnt from the ITPA-IOS joint activity of plasma burn-through simulation code benchmark [17], we have developed full circuit equations and plasma volume model in DYON. The mathematical representation of the improved model has been provided in this paper. The improved model enables plasma burn-through simulations only with the control room operation parameters such as coil currents (or voltages) in solenoid and PF coils, prefill gas pressure, and gas fuelling data. This new feature is necessary to perform predictive simulations and optimize operation scenarios in a future device such as ITER and STEP. In addition, the plasma volume evolution during the burn-through phase can now be consistently calculated with the physics-based model, reducing one of the main uncertainties that plasma burn-through modelling has had so far. The upgraded plasma burn-through model has successfully reproduced the 34 MAST discharge with the direct induction scenario, indicating the validity of the upgraded model and the prediction capability. For more complete validation, dedicated experiments with operation parameter scanning and measurement of plasma volume with open field configuration are needed. It will improve the confidence when using the electromagnetic plasma burn-through model to optimize ITER operation scenarios and to design STEP.