Metis: a fast integrated tokamak modelling tool for scenario design

Among the quantities simulated by core transport codes, the plasma current density is the one that is most accurately predicted in a large variety of experiments. In the absence of magneto-hydrodynamic activity, the neoclassical resistivity seems a valid diffusive model to describe resistive current diffusion in tokamak experiments, at least during the flat-top phase of the H-mode [9]. Therefore in METIS current diffusion and plasma equilibrium computations are fully kept in the 1.5D paradigm: (i) a 2D equilibrium code is run periodically to update the metrics of the poloidal flux transport equation in a way consistent with the plasma pressure and current profiles; (ii) the 1D poloidal flux transport equation is solved with all its terms.

The current diffusion equation is solved in terms of the poloidal flux Ψ, on a uniform toroidal flux coordinate ρ grid, exactly as in the CRONOS 1.5D code [2]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\left. \frac{\partial \Psi }{\partial t} \right|}_{{{\rho }_{{\rm norm}}}}}=&\frac{\left\langle \frac{{{\left| \nabla \rho \right|}^{2}}}{{{R}^{2}}} \right\rangle }{{{\mu }_{0}}{{\sigma }_{\parallel }}\rho _{{\rm m}}^{2}\left\langle \frac{1}{{{R}^{2}}} \right\rangle }\frac{{{\partial }^{2}}\Psi }{\partial \rho _{{\rm norm}}^{2}}+\left. \Bigg\{\frac{\left\langle \frac{{{\left| \nabla \rho \right|}^{2}}}{{{R}^{2}}} \right\rangle }{{{\mu }_{0}}{{\sigma }_{\parallel }}\rho _{{\rm m}}^{2}\left\langle \frac{1}{{{R}^{2}}} \right\rangle }\frac{\partial }{\partial {{\rho }_{{\rm norm}}}}\left[ \ln \left( \frac{\frac{\partial V}{\partial {{\rho }_{{\rm norm}}}}\left\langle \frac{{{\left| \nabla \rho \right|}^{2}}}{{{R}^{2}}} \right\rangle }{F} \right) \right]\right.\nonumber \\ &\left.+~\frac{{{\rho }_{{\rm norm}}}}{{{\rho }_{{\rm m}}}}\frac{{\rm d}{{\rho }_{{\rm m}}}}{{\rm d}t}+\frac{{{\rho }_{{\rm norm}}}}{2{{B}_{0}}}\frac{{\rm d}{{B}_{0}}}{{\rm d}t} \Bigg.\right\}\frac{\partial \Psi }{\partial {{\rho }_{{\rm norm}}}}+~\frac{{{B}_{0}}}{{{\sigma }_{\parallel }}F\left\langle \frac{1}{{{R}^{2}}} \right\rangle }~{{j}_{{\rm ni}}}\nonumber \end{align} \tag{ 1 }$$

where ${{\left. \frac{\partial \Psi }{\partial t} \right|}_{{{\rho }_{{\rm norm}}}}}$ is the time derivative of the poloidal flux at a given radial position ${{\rho }_{{\rm norm}}}$, σ_|| denotes the parallel conductivity (calculated following the Sauter model [10]), F is the diamagnetic function, j_ni the current density driven by the non-inductive sources, R the major radius, ${{\mu }_{0}}$ the magnetic permeability of free space, (4π  ×  10⁻⁷ in MKS units), ${{\rho }_{{\rm m}}}$ the value of ρ at the last closed flux surface, and the normalised toroidal flux coordinate ${{\rho }_{{\rm norm}}}=\frac{\rho }{{{\rho }_{{\rm m}}}}$. The notation 〈〉 indicates a magnetic flux surface average, defined as the volume average of a quantity around a flux surface of radial coordinate ρ, i.e. in an elementary volume dV enclosed between two magnetic surfaces distant of dρ. We remind the definition of the toroidal flux coordinate, $\rho =\sqrt{\frac{\Phi }{\pi {{B}_{0}}}}$, where Φ is the toroidal magnetic flux, B₀ the vacuum magnetic field at a given major radius R₀ (usually taken at the centre of the vacuum vessel). The normalised radial grid does not depend on time, but ρ_m is time-dependent and calculated by the equilibrium solver.

Solving the current diffusion in the same way as the traditional 1.5D integrated modelling code is not a drawback from the performance point of view since it corresponds to the transport phenomenon with the longest time scale. To speed up the calculation, current diffusion is calculated on a 21 points radial grid only and the 3-moments description (Shafranov Shift, ellipticity and triangularity) is used for the MHD equilibrium [11, 12]. The last closed flux surface (LCFS) can be either described by moments or as a series of (R, Z) points. In this last case an additional morphing (continuous deformation) is applied to flux surfaces to match the LCFS with effect weighted as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho \left(x,\theta \right)=~{{\rho }_{{\rm moment}}}\left(x,\theta \right)~\left[ \left(1-~{{x}^{{{e}_{{\rm m}}}}} \right)+~{{x}^{{{e}_{{\rm m}}}}}~\frac{{{\rho }_{{\rm LCFS}}}\left(\theta \right)}{{{\rho }_{{\rm moment}}}\left(1,\theta \right)} \right] \nonumber \end{align} \tag{ 2 }$$

where ${{\rho }_{{\rm moment}}}\left(x,\theta \right)$ is the moment description of flux surfaces ($x$ is the radial coordinate described in next section and $\theta$ is a poloidal angle); ${{\rho }_{{\rm LCFS}}}\left(\theta \right)$ is the LCFS poloidal representation given as a series of (R, Z) points; ${{\rho }_{{\rm moment}}}\left(1,\theta \right)$ is the LCFS described by moment and ${{e}_{{\rm m}}}$ is an arbitrary parameter. The value ${{e}_{{\rm m}}}=5$ is typically used, as it has been found to minimize the discrepancy for geometrical coefficients involved in the current diffusion equation (1) for various devices (JET, Tore Supra, JT-60SA, ITER) and scenarios, with respect to equilibria calculated using the HELENA code [13]. The value of this parameter is kept tunable by the user.

For the description of all radial profiles, METIS uses a uniform 21 points grid on the normalized minor radius x, defined as follows: for each magnetic surface, the minor radius is defined as $a=\frac{R\max -R\min }{2}$, where R_max (resp. R_min) is the maximum (resp. minimum) major radius of that surface. The minor radius is then normalized to its value at the LCFS a_m, so that $x=\frac{a}{{{a}_{{\rm m}}}}$.

In H-mode, the top of pedestal, located at ${{x}_{{\rm ped}}}$ is always set in METIS at the second point from the edge. The internal grid of METIS is ${{x}_{k}}=\frac{{{a}_{{\rm m}}}}{20}\left(k-1 \right)$ with $k$ an integer between 1 and 21 and ${{x}_{{\rm ped}}}=~{{x}_{20}}$. Therefore, the width of the pedestal is always 5% of plasma minor radius. Use of a simple grid with fixed step simplifies the numerical computation. This technical choice typically overestimates the pedestal width which is about 2%–3% of the minor radius in present experiments, but this has a very limited impact on the profiles predicted by METIS, which are based on scaling expressions for the energy content of core plasma and the pedestal (so they are not gradient based, see the section below). We have checked that with a width of 5% we have about the same amount of bootstrap current (the difference is less than 5%), integrated on pedestal width in METIS than with a pedestal width of 2%–3% in CRONOS simulations.

The transport of energy in the core of tokamak plasmas is dominated by turbulence. In spite of enormous efforts and progress made by the fusion community to understand and predict plasma turbulence, the prediction of the resulting transport flux still features large uncertainties. First-principle gyro-kinetic codes are extremely expensive in terms of computing resources so that they can barely be applied to the long-time scales of plasma experiments. On the other hand, the existing reduced models (e.g. GLF23 [14], TGLF [15], Qualikiz [16]), even when first-principle based, fail to capture completely the complexity of the turbulence phenomena and their predictions are in some cases quite far from the experimental results. This is particularly true for plasmas at high beta and a significant fraction of energetic particles, where non-linear physics becomes dominant [17, 18], but also for H-mode plasmas where the fusion performance is fully dominated by the height of the pedestal which is still far from being predicted reliably. We touch here the main paradox of present Integrated Modelling, which is that the largest source of uncertainties lies at the heart of the problem solved, namely the prediction of the heat and particle fluxes. Moreover, the reduced 'first-principle based' models usually feature strong dependences on the gradients of the transported quantities, which requires using rather small time steps in the transport solver (typically ~10⁻⁴ s) if one wants a numerically accurate resolution of the transport dynamics. This numerical constraint makes this type of models difficult to apply to simulate a full ITER discharge lasting several hundreds of seconds. This is even more difficult in view of scenario optimisation studies where it is desired to test several combinations of actuators, e.g. power, timing and parameters, hence requiring a large number of these long simulations and some manual trial and error process.

Therefore simplifying the heat transport modelling assumptions, which are both the largest source of uncertainty and the largest CPU time consumption in a classical 1.5D integrated modelling simulation, is an effective way of increasing the performance with the smallest impact on the reliability of the results, in view of a fast scenario simulator. In METIS, heat transport is treated in a mixed 0D–1D approach in two steps, by separating the time and radial dimensions. This approach is found to be quite successful for simulating the main dynamics of the core plasma temperature profiles, although it naturally has limitations for phenomena in which spatial and temporal evolutions are coupled: for example, details of heat pulse propagation following a pellet injection cannot be resolved and would require using a classical 1.5D solver.

The first step consists in solving a time-dependent 0D equation for the plasma thermal energy content W_th:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{{\rm d}{{W}_{{\rm th}}}}{{\rm d}t}=-\frac{{{W}_{{\rm th}}}}{{{\tau }_{{\rm E}}}}+{{P}_{{\rm loss}}}.\nonumber \end{align} \tag{ 3 }$$

In this ordinary differential equation scaling expressions for the energy confinement time τ_E are typically used, a multi-machine approach widely used for extrapolating present results to future tokamaks such as ITER [19]. Although this approach has inherent limitations when non-linear physics is dominant, it is still valid for a significant class of future plasmas [20]. P_loss represents the total power transported through the plasma separatrix by diffusion or convection mechanisms (detailed expressions are given in appendix B). Multiple scaling expressions are coded in METIS for calculating the value of τ_E in a self-consistent way with the actual parameters. The expression of W_th actually depends on whether the plasma is in L or H mode. This information, as well as the ability to handle L–H and H–L transitions, is of course essential for a full scenario simulator such as METIS. The list of available scaling expressions is of course easily extensible and could cope with e.g. new scaling laws based on the future analysis of the ITER shot database, in order to make METIS even more useful for ITER operation.

The METIS code uses scaling expressions of the L–H power threshold to deduce whether the plasma is in L or H mode. Multiple options are available from the literature and new ones can be easily added.

When in L mode, no pedestal is created and an L-mode scaling expression is used for τ_E in equation (3).

When in H mode:

• A pedestal is created. Its height can be prescribed in multiple ways (constant or scaling expression).
• An H-mode scaling expression is used for τ_E in equation (1).

By default, the L–H transition is modeled as an immediate change of the value of τ_E in equation (1). Nevertheless in experiments, the transition from L to H mode is not abrupt when crossing the power threshold. There is often a P_loss range slightly above the threshold P_threshold that yields an intermediate confinement level. This can be mimicked in METIS by defining a linear transition in $\frac{{{P}_{{\rm loss}}}}{{{P}_{{\rm threshold}}}}$ between the L and H mode energy confinement scaling expressions. A detailed explanation about how P_loss is calculated is given in the appendix.

The second step consists in calculating the electron and ion temperature profiles, assuming steady-state transport equations and purely diffusive transport, the heat flux at a radial position x can be expressed as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial {{T}_{{\rm e}}}}{\partial x}=\frac{-\int_{0}^{x}{{{V}^{\prime }}{{Q}_{e}}}}{{{n}_{{\rm e}}}{{\chi }_{{\rm e}}}{{V}^{\prime }}\langle \mid \nabla \rho {{\mid }^{2}}\rangle }\,{\rm and}\,\frac{\partial {{T}_{{\rm i}}}}{\partial x}=\frac{-\int_{0}^{x}{{{V}^{\prime }}{{Q}_{i}}}}{{{n}_{{\rm i}}}{{\chi }_{{\rm i}}}{{V}^{\prime }}\langle \mid \nabla \rho {{\mid }^{2}}\rangle }.\nonumber \end{align} \tag{ 4a }$$

Or alternatively the conductivity can be used instead of diffusivity:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial {{T}_{{\rm e}}}}{\partial x}=\frac{-\int_{0}^{x}{{{V}^{\prime }}{{Q}_{{\rm e}}}}}{{{\kappa }_{{\rm e}}}{{V}^{\prime }}\langle \mid \nabla \rho {{\mid }^{2}}\rangle }\,{\rm and}\,\frac{\partial {{T}_{{\rm i}}}}{\partial x}=\frac{-\int_{0}^{x}{{{V}^{\prime }}{{Q}_{{\rm i}}}}}{{{\kappa }_{{\rm i}}}{{V}^{\prime }}\langle \mid \nabla \rho {{\mid }^{2}}\rangle }\nonumber \end{align} \tag{ 4b }$$

where Q_e and Q_i are the sum of all the electron and ion heat source terms, including the equipartition term Q_ei. The diffusion coefficients for electron and ions are noted respectively ${{\chi }_{{\rm e}}}$ and ${{\chi }_{{\rm i}}}$. The conductivity coefficients for electron and ions are noted respectively ${{\kappa }_{{\rm e}}}$ and ${{\kappa }_{{\rm i}}}$ (linked to diffusion coefficients by: ${{\kappa }_{{\rm e},{\rm i}}}={{n}_{{\rm e},{\rm i}}}{{\chi }_{{\rm e},{\rm i}}}$). The geometrical coefficient ${{V}^{\prime }}$ is the derivative of the plasma volume enclosed in a magnetic surface with respect to the normalised minor radius x, while $\langle \mid \nabla \rho {{\mid }^{2}}\rangle$ is the surface average of the squared gradient of the toroidal flux coordinate. We note that although the temperature profiles are solved using steady-state equations, the dynamics of heat transport can still be accounted for by the first equation (on the global energy) which is time-dependent and used to normalize the heat conductivities ${{\kappa }_{{\rm e}}}$ and ${{\kappa }_{{\rm i}}}$ (or alternatively the diffusion coefficients ${{\chi }_{{\rm e}}}$ and ${{\chi }_{{\rm i}}}$). Validation against experiments shows that this approximated approach is relevant for the description of transient phenomena at the scenario level (e.g. transport dynamics during current ramps). Note also that with this formulation, the temperature profiles calculated by METIS take into account the information on the radial distribution of the heat sources.

To complete the second step, the diffusion coefficients used in the equations (4) must be calculated in a way consistent with the energy content that has been calculated at the first step. This is done differently in L and H mode.

In L-mode:

• i.  
  The dependence of the electron diffusion coefficient on plasma parameters or on the radial coordinate is prescribed. Multiple options are coded in METIS and new expressions can be easily added. Three widely used options are (a) Bohm-gyro-Bohm model [21, 22] (b) fixed radial dependence of the type ${{\chi }_{{\rm e}}}={{\chi }_{0}}(1+\lambda {{x}^{\gamma }})$; where $\lambda$ and $\gamma$ are constants chosen by the user (the default values are $\lambda =3$ and $\gamma =2$, or (c)${{\chi }_{{\rm e}}}={{\chi }_{0}}{{q}^{\gamma }}\left(x \right)$; where q is the safety factor and $\gamma$ is a constant chosen by the user.
• ii.  
  The ion diffusion coefficient is calculated from the electron one with a simple expression of the type: ${{\chi }_{{\rm i}}}={{\mu }_{{\rm e},{\rm i}}}{{\chi }_{{\rm e}}}$. Often ${{\mu }_{{\rm e},{\rm i}}}$ is prescribed to be a constant, but it can also be provided by other analytical expressions [23–26]
• iii.  
  Steps i and ii yield the radial shape of the diffusion coefficients and their relative value. The coefficients are then normalized in such a way that the thermal energy content obtained from the integral of the pressure profile corresponds to the value calculated from equation (1), i.e. $\frac{3}{2}\int_{x=0}^{1}{({{n}_{{\rm e}}}{{T}_{{\rm e}}}+{{n}_{{\rm i}}}{{T}_{{\rm i}}}){{V}^{\prime }}{\rm d}x={{W}_{{\rm th}}}}$. In order to have an exact (inverse) proportionality relation between the normalization factor of the diffusivities involved in equations (4) and the energy content, the contribution of the boundary conditions ${{T}_{{\rm e},{\rm LCFS}}}$ and ${{T}_{{\rm i},{\rm LCFS}}}$ to the integral must be removed. Indeed equations (4) define only the temperature gradient, therefore the temperature profiles are defined as $T(x)={{T}_{{\rm LCFS}}}+\int_{1}^{x}{\frac{\partial T}{\partial x}{\rm d}x}$, in which only the second term is inversely proportional to the diffusivity. The contribution to the integral of a finite LCFS temperature represents an energy content ${{W}_{0}}=\frac{3}{2}$$\int_{x=0}^{1}{({{n}_{{\rm e}}}{{T}_{{\rm e},{\rm LCFS}}}+{{n}_{{\rm i}}}{{T}_{{\rm i},{\rm LCFS}}}){{V}^{\prime }}{\rm d}x}$. The normalization constant of the diffusivities is therefore calculated following the relation $\frac{3}{2}\int_{x=0}^{1}{({{n}_{{\rm e}}}\int_{1}^{x}{\frac{\partial {{T}_{{\rm e}}}}{\partial x}{\rm d}x+{{n}_{i}}}}\int_{1}^{x}{\frac{\partial {{T}_{{\rm i}}}}{\partial x}{\rm d}x){{V}^{\prime }}{\rm d}x={{W}_{{\rm th}}}-{{W}_{0}}}$.

To complete step 2, a convergence loop on the temperature profiles resulting from this process is performed in order to find the self-consistent value of the equipartition term Q_ei in equations (4), which is proportional to T_e  −  T_i.

In H-mode, the procedure is similar to that of the L mode, with the only difference that the pedestal top is used as the boundary condition for the integration of equations (4). A scheme of such a procedure is shown in figure 1.

Zoom In Zoom Out Reset image size

Particle transport is also dominated by turbulence, however the level of understanding is lower than that for heat transport. The existence of significant inward and outward pinches which in turn depend on the accurate description of turbulence requires a significant amount of computational time if a first principle approach is used. Additionally, unlike for heat transport, the sources uncertainties are significant close to the pedestal region. On the other hand, a significant amount of work has been performed in order to characterize some density features, as the peaking, by analysing extensive plasma tokamak databases [27]. In METIS, an approach that is even simpler than for heat transport is chosen, based on the prescription of some quantities.

The primary quantity for determining the density of the plasma species in METIS is the electron density profile. It is described by 3 parameters, which can vary with time:

• i.  
  the line-averaged density (${{\bar{n}}_{{\rm e}}}$), which is prescribed
• ii.  
  the peaking factor $\frac{{{n}_{0}}}{\left\langle n \right\rangle }$ which is either prescribed or computed with the help of a scaling expression (to take into account self-consistent dependencies on plasma parameters, various options from the literature are available, with one set of scaling expression for L-mode and another one for H-mode)
• iii.  
  the density value at the separatrix (${{n}_{{\rm e},a}}$) obtained from simple models or scaling expressions depending of the plasma configuration: poloidal limiter, toroidal limiter, axisymmetric divertor with X-point. In the last case (X-point) the expression depends also of the confinement mode L or H.

In L mode the shape of the profile is defined as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{n}_{{\rm e}}}(t,x)=({{n}_{{\rm e},0}}(t)-{{n}_{{\rm e},a}}(t)){{(1-{{x}^{2}})}^{{{\nu }_{n}}(t)}}+{{n}_{{\rm e},a}}(t)\nonumber \end{align} \tag{ 5 }$$

where ${{\nu }_{n}}=\frac{{{n}_{{\rm e}}}(t,x=0)}{\langle {{n}_{{\rm e}}}\rangle (t)}-1$.

In H mode, in order to ensure that the pedestal is also present in the electron density profile, the density profile is computed by another method. The density profile is constrained by the line-averaged value, the peaking factor, the constraint of flat profile at the centre, the edge value and the constraint that the temperature profiles are strictly monotonically decreasing. A piecewise cubic Hermite polynomial interpolation is used to compute the profile, calculating the pedestal density that fulfils those 4 constraints: the electron density is given at 3 points: at the centre (x  =  0), at the top of the pedestal (x  =  0.95) and at the edge (x  =  1) with the constraint of null derivative at the centre. The value at the top of the pedestal is computed as the maximal value ensuring negative derivative of electron temperature (with a minimum of $0.01~\left({{T}_{{\rm e},0}}-{{T}_{{\rm e},a}} \right)$)

From the results of an initial METIS simulation, it is possible to compute as a post-processing the electron temperature, ion temperature and electron density using more sophisticated transport models than the scaling-based ones. For this post-processing duration to remain of the same order as the one of the main METIS computation, neural network based models such as Qualikiz-NN should be used, which are approximations to quite sophisticated transport models [28]. This allows comparing the prediction of such more physics-based transport models to the scaling based-ones. In this calculation, the steady-state transport equation is solved from the top of the pedestal to the magnetic axis and keeping the sources, current density profile and equilibrium from the result of the initial METIS simulation. The steady-state solution is thus calculated for each time step using physics-based transport models. However, in order to keep a fast computation time, this is at the expense of losing the self-consistency between transport, equilibrium, and sources.

In METIS, all ion species density profiles are assumed to be proportional to the electron density profile, with exceptions of the tungsten species and helium ashes, which have specific treatments (see appendices D.1 and D.2). Alternatively to this simple assumption, one can use a simple neoclassical model based on impurity shielding/accumulation, similar to that one used for tungsten. However in practice this model is not frequently used, as turbulent transport remains generally dominant over the neoclassical one for light impurities.

The user also specifies the reference effective charge (line averaged) for the plasma Z_eff,ref, which is either prescribed as a time-dependent value or self-consistently calculated from a scaling expression. This value is called 'reference', although it is used as the exact effective charge of the plasma when the main plasma species are H, D, He. A correction is applied in the case of D–T mixture to account for the presence of ${{H}_{{\rm e},{\rm ashes}}}$ resulting from the D–T fusion reactions and in case of presence of tungsten impurity. It thus represents the effective charge in absence of ${{H}_{{\rm e},{\rm ashes}}}$ and tungsten.

Then, various types of plasma compositions, which have to be consistent with Z_eff,ref, can be prescribed by the user. METIS first calculates their volume average densities, from the rules explained below, then applies the radial profile option (by default: same profile as the electron density) to these average densities.

METIS distinguishes the following types of ion species:

• Main species: can be H, D, D–T mixture, or He. In case of a D–T mixture, the user prescribes the isotopic density ratio $\frac{\left\langle {{n}_{{\rm T}}} \right\rangle }{\left\langle {{n}_{{\rm D}}} \right\rangle }$ and a specific treatment is applied for the Helium density (see section D.1).
• Minority species: in case of ion cyclotron resonant heating (ICRH) minority heating, the user can specify the type of one minority species (H, D,T, He3 and He4) and the average density ratio of the minority species with respect to the main ion species.
• Impurity species: two additional impurity species can be included, the user specifies their type as well as their relative average density ratio $\frac{\left\langle {{n}_{{\rm imp},1}} \right\rangle }{\left\langle {{n}_{{\rm imp},2}} \right\rangle }$.
• Optionally, the tungsten species can be added. It is treated separately, self consistently with the divertor to provide feedback on the tungsten source (see section D.2).

2.6.1. Computation of the ion densities.

The specification of the plasma composition and the average density ratios as described above, together with the effective charge and electro-neutrality constraints, provides a linear system with equal number equations and unknowns and thus allows calculating the average density of all ion species.

The final value of the effective charge is computed with the whole set of profiles including helium ashes and tungsten.

In the presence of Helium ashes or Tungsten impurity, the effective charge is modified and become:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{Z}_{{\rm eff}}}\left(t,x \right)=~\frac{{{n}_{{\rm H}}}\left(t,x \right)+{{n}_{{\rm D}}}\left(t,x \right)+{{n}_{{\rm T}}}\left(t,x \right)+4~{{n}_{{\rm He}}}\left(t,x \right)+~Z_{{\rm imp},1}^{2}{{n}_{{\rm imp},1}}\left(t,x \right)+Z_{{\rm imp,}2}^{2}{{n}_{{\rm imp},2}}\left(t,x \right)~+~Z_{{\rm W},{\rm averaged}}^{2}\left({{T}_{{\rm e}}}\left(t,x \right) \right){{n}_{{\rm W}}}\left(t,x \right)}{{{n}_{{\rm e}}}\left(t,x \right)}.\nonumber \end{align} \tag{ 6 }$$

Verifying (in case of D–T plasma):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{n}_{{\rm e}}}\left(t,x \right)=&~{{n}_{{\rm H}}}\left(t,x \right)+{{n}_{{\rm D}}}\left(t,x \right)+{{n}_{{\rm T}}}\left(t,x \right)\nonumber \\ &+4~{{n}_{{\rm He}}}\left(t,x \right)+~{{Z}_{{\rm imp},1}}{{n}_{{\rm imp},1}}\left(t,x \right)\nonumber \\ &+{{Z}_{{\rm imp},2}}~{{n}_{{\rm imp},2}}\left(t,x \right)+{{Z}_{{\rm W},{\rm averaged}}}\left({{T}_{{\rm e}}}\left(t,x \right) \right){{n}_{{\rm W}}}\left(t,x \right).\nonumber \end{align}$$

With ${{n}_{{\rm imp},2~=~}}\frac{\left\langle {{n}_{{\rm imp},1}} \right\rangle }{\left\langle {{n}_{{\rm imp},2}} \right\rangle }{{n}_{{\rm imp},1}}$, ${{n}_{{\rm T}}}=\frac{\left\langle {{n}_{{\rm T}}} \right\rangle }{\left\langle {{n}_{{\rm D}}} \right\rangle }~{{n}_{{\rm D}}}$ and ${{n}_{{\rm H}}}={{c}_{\min }}{{n}_{{\rm D}}}$

For other main ions choice the formulation is updated as follows:

• Hydrogen main ion: if the minority species is not deuterium ${{n}_{{\rm D}}}=\frac{\left\langle {{n}_{{\rm T}}} \right\rangle }{\left\langle {{n}_{{\rm D}}} \right\rangle }{{n}_{{\rm H}}}$ and ${{n}_{{\rm T}}}=0$, otherwise ${{n}_{{\rm D}}}=\bigg (\frac{\left\langle {{n}_{{\rm T}}} \right\rangle }{\left\langle {{n}_{{\rm D}}} \right\rangle }+$${{c}_{\min }} \bigg)~{{n}_{{\rm H}}}$ and ${{n}_{T}}=~0$.
• Helium main ion: ${{n}_{{\rm D}}}={{n}_{{\rm H}}}={{n}_{{\rm T}}}=0$; if ICRH minority heating scheme using H, D or T then ${{n}_{{\rm H},{\rm D},{\rm T}}}={{c}_{\min }}{{n}_{{\rm He}}}$

An estimate of the rotation is carried out in METIS, mainly taking into account the effect of neutral beam injection and intrinsic rotation. Simple models are used to handle toroidal rotation due to parallel electric field and lower-hybrid current drive (LHCD). The effect of ripple losses is only taken into account for Tore Supra, device for which a scaling exists [29]. The toroidal rotation due to fast ion losses and to fast ion momentum transport are not taken into account, because there is no simple model available for these. No feedback of the rotation on the standard confinement is taken into account. Toroidal rotation only impacts the ITB formation.

The simple model implemented in METIS is sufficient for the study of neutral beam injection (NBI) dominated plasmas, as well as for reactor studies, in which the fast alpha distribution is close to be isotropic and does not allow to transport a significant part of toroidal momentum [30]. This simplified model does not allow studying plasma rotation when the plasma is heated mainly by ion cyclotron resonant heating (ICRH), when ripple is significant on other devices than Tore Supra, or when the confinement of fast particles is poor (i.e. when the fast particle orbit size is comparable to the minor radius of the plasma).

The computation of rotation in METIS is separated in three parts. The first part consists in the evaluation of the volume-averaged toroidal momentum. The second part consists in the evaluation of rotation at the LCFS. The third part consists in the computation of the radial profiles of toroidal and poloidal rotation velocities and finally of the radial electric field.

The volume-averaged toroidal momentum computation is analogous to that of the energy content. A general discussion can be found in [31]. The total momentum is defined as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{R}_{{\rm tot}}}=\int_{0}^{1}{\underset{k\in \left\{{\rm species} \right\}}{\mathop \sum }\,{{m}_{{\rm p}}}{{A}_{k}}{{n}_{k}}R{{v}_{\varphi ,k}}{{V}^{\prime }}{\rm d}x}\nonumber \end{align} \tag{ 7 }$$

where ${{v}_{\varphi ,k}}$ is the toroidal velocity of the species k, ${{m}_{{\rm p}}}$ the proton mass, ${{A}_{k}}$ the number of nucleons of the species, ${{n}_{k}}$ the species density, $R$ the geometrical major radius of the flux surface. With this approach, the contribution of electrons is neglected because of their negligibly small mass.

The evolution equation of ${{R}_{{\rm tot}}}$ is:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{{\rm d}{{R}_{{\rm tot}}}}{{\rm d}t}=-\frac{{{R}_{{\rm tot}}}}{{{\tau }_{\varphi }}}+{{S}_{\varphi ,{\rm NBI}}}+{{S}_{\varphi ,{\rm self}}}+~{{S}_{\varphi ,RF}}+{{S}_{\varphi ,{{E}_{\parallel }}}}+~{{S}_{\varphi ,{\rm ripple}}}+{{F}_{n,0}}\nonumber \end{align} \tag{ 8 }$$

where the toroidal momentum confinement time ${{\tau }_{\varphi }}$ is assumed to be related to the energy confinement time ${{\tau }_{{\rm E}}}$ (see below for details) and includes the friction of the plasma with neutrals at the edge. From ${{R}_{{\rm tot}}}$ we deduce the volume-averaged plasma rotation ${{\omega }_{\varphi }}=\frac{{{R}_{{\rm tot}}}}{{{\Gamma }_{\varphi }}{{R}_{{\rm axis}}}}$ with ${{\Gamma }_{\varphi }}$ is the conversion factor between velocity and momentum:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\Gamma }_{\varphi }}=\frac{1}{\left\langle {{v}_{\varphi ,{\rm shape}}} \right\rangle }\int_{0}^{1}{\underset{k\in \left\{{\rm species} \right\}}{\mathop \sum }\,{{m}_{{\rm p}}}{{A}_{k}}{{n}_{k}}R{{v}_{\varphi ,{\rm shape}}}{{V}^{\prime }}{\rm d}x}\nonumber \end{align} \tag{ 9 }$$

and $\left\langle {{v}_{\varphi ,{\rm shape}}} \right\rangle =~\frac{\int_{0}^{1}{{{v}_{\varphi ,{\rm shape}}}{{V}^{\prime }}}{\rm d}x}{\int_{0}^{1}{{{V}^{\prime }}}{\rm d}x}$

with ${{v}_{\varphi ,{\rm shape}}}$ is the profile of toroidal rotation as described below (assuming that all species have same rotation shape profile) and R_axis the center of each flux surface.

Momentum sources are due to:

• 1.  
  Neutral beam injection toroidal moment source:

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{S}_{\varphi ,{\rm NBI}}}=\sum\limits_{{\rm b}\in \left\{{\rm H},{\rm D},{\rm T} \right\}}{{{m}_{{\rm b}}}\int_{0}^{1}{{{R}_{{\rm axe}}}\left(x \right)\frac{{{p}_{{\rm NBI},{\rm b}}}\left(x \right)}{e{{E}_{{\rm b}0}}}\sqrt{\frac{2e{{E}_{{\rm b}0}}}{{{m}_{{\rm b}}}}}~{{\mu }_{{\rm b}}}\left(x \right){{V}^{\prime }}\left(x \right)}}{\rm d}x.\nonumber \end{align} \tag{ 10 }$$

  With ${{R}_{{\rm axe}}}$ is the magnetic axis of each flux surface; ${{p}_{{\rm NBI},{\rm b}}}\left(x \right)$ is the profile of power deposition due to NBI; ${{\mu }_{{\rm b}}}\left(x \right)$ is the profile of pitch angle ($\mu =~\frac{{{v}_{\parallel }}}{v}$); ${{m}_{{\rm b}}}$ is the mass of injected species; ${{E}_{{\rm b}0}}$ is the injection energy in eV and e the electron charge.This expression implies momentum conservation, which is true if fast ion losses are negligible.
• 2.  
  The intrinsic plasma toroidal rotation source is:

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{S}_{\varphi ,{\rm self}}}=\frac{{{\Gamma }_{\varphi ~{{v}_{\varphi ,{\rm self}}}}}}{{{\tau }_{\varphi }}}.\nonumber \end{align} \tag{ 11 }$$

  The intrinsic plasma rotation ${{v}_{\varphi ,{\rm self}}}$ is given by a scaling law time a factor:

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{v}_{\varphi ,{\rm self}}}=~{{f}_{{\rm intrinsic}}}~{{v}_{\varphi ,{\rm scaling}}}\nonumber \end{align} \tag{ 12 }$$

  Where ${{v}_{\varphi ,{\rm scaling}}}$ is either:

  • Rice scaling is taken from [32]:
    $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{v}_{\varphi ,{\rm scaling}}}={{0.4610}^{11}}B_{{\rm ref}}^{1.1}P_{\varphi }^{1.0}I_{{\rm p}}^{-1.9}R_{{\rm ref}}^{2.2}\nonumber \end{align} \tag{ 13 }$$
    $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{P}_{\varphi }}=2\left\langle e{{T}_{{\rm i}}}{{\sum }_{k\in \left\{{\rm species} \right\}}}{{n}_{k}} \right\rangle +\frac{2}{3}\frac{\left(W-{{W}_{{\rm th}}} \right)}{{{V}_{{\rm p}}}}.\nonumber \end{align} \tag{ 14 }$$
    This expression for ${{P}_{\varphi }}$ takes into account the fact that for dominant electron heating the spontaneous toroidal rotation changes only because of ion pressure increased. On JET it was found that for NBI heated plasma, the thermal pressure is close to two times the ion pressure.
  • Barnes and Parra scaling [33] for which we have taken $\Delta {{T}_{{\rm i}}}=~\frac{\left\langle {{P}_{{\rm ion}}} \right\rangle }{\left\langle {{n}_{{\rm ion}}} \right\rangle }$
  • deGrassie scaling [34] for which we have taken $\Delta {{T}_{{\rm i}}}=~\frac{\left\langle {{P}_{{\rm ion}}} \right\rangle }{\left\langle {{n}_{{\rm ion}}} \right\rangle }$.

  This definition of $\Delta {{T}_{{\rm i}}}$ is the best generalization of the experimental measurement used in papers that can extend the implementation of the formula to METIS.The factor ${{f}_{{\rm intrinsic}}}$ is either provided by the user or computed using the model described in [35].
• 3.  
  RF driven rotation source terms:This source term correspond to the toroidal moment transferred from RF waves to the plasma [36]:

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{S}_{\varphi ,{\rm RF}}}=~-~(2~{{D}_{{\rm LH}}}-1)~\frac{{{R}_{{\rm ref}}}~{{n}_{\parallel ,0}}}{c}~{{P}_{{\rm LH},{\rm th}}}-{{s}_{{\rm ECCD}}}~\frac{{{R}_{{\rm ref}}}}{2~c}~{{P}_{{\rm ECRH}}} \nonumber \end{align} \tag{ 15 }$$

  where ${{D}_{{\rm LH}}}$ is the directivity of LH wave. The parallel refractive index of the EC wave is assumed to be near the optimal value to maximized current drive efficiency (${{s}_{{\rm eccd}}}\sim 0.5$) ${{R}_{{\rm ref}}}$ is the geometrical center of the plasma (${{R}_{{\rm ref}}}=\frac{{{R}_{\max ~}}-{{R}_{\min }}}{2}$).
• 4.  
  Rotation source due to parallel electric field:This term is negligible during flat-top but important during ramp-up and also because it breaks the symmetry between co and counter current [37]. This term reads:

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{S}_{\varphi ,{{E}_{\parallel }}}}=\frac{{{m}_{{\rm e}}}}{e~{{m}_{{\rm p}}}}\frac{{{\Gamma }_{\varphi }}}{{{\tau }_{\varphi }}}\frac{{{M}_{{\rm eff}}}}{{{Z}_{{\rm eff}}}}\frac{{{I}_{\Omega }}+{{I}_{{\rm run}}}}{\left\langle {{n}_{e}}~ \right\rangle {{S}_{{\rm p}}}} \nonumber \end{align} \tag{ 16 }$$

  where ${{M}_{{\rm eff}}}$ is the effective mass of the plasma; ${{I}_{\Omega }}$ the ohmic current and ${{I}_{{\rm run}}}$ the electron runaway current; $\left\langle {{n}_{{\rm e}}} \right\rangle$ is the volume averaged electron density and ${{S}_{{\rm p}}}$ the plasma poloidal surface.
• 5.  
  Friction on cold neutral:At the edge of the core plasma, the rotation is slowed down by the friction with cold neutrals coming from recycling and gas puff. The main contribution comes from charge exchange between cold neutrals and plasma ions. This term reads:

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{F}_{n0}}=-{{m}_{{\rm p}}}\int_{0}^{1}{\left(\sum\limits_{k\in \left\{{\rm H},{\rm D},{\rm T} \right\}}{{{A}_{k~}}{{n}_{k}}}~ \right){{\left\langle \sigma v \right\rangle }_{{\rm cx}}}{{R}_{{\rm axe}}}{{v}_{\varphi ~}}{{n}_{0}}{{V}^{\prime }}}{\rm d}x\nonumber \end{align} \tag{ 17 }$$

  where ${{n}_{0}}$ is the density of neutral hydrogen isotopes (see appendix A.11), ${{\left\langle \sigma v \right\rangle }_{{\rm cx}}}$ the rate for charge exchange reactions. This factor is implemented only for H/D/T plasma.

The confinement time of toroidal rotation is observed to be a fraction of the energy confinement time [38, 39] and often close to the ion confinement time:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\tau }_{\varphi }}={{f}_{\tau ,{\rm rot}}}\min \left({{\tau }_{{\rm E}}},{{\tau }_{{\rm ii}}} \right).\nonumber \end{align} \tag{ 18 }$$

With ${{\tau }_{{\rm ii}}}=\frac{e\int_{0}^{1}{{{n}_{{\rm i}}}{{T}_{{\rm i}}}{{V}^{\prime }}{\rm d}x}}{\int_{0}^{1}{\left({{Q}_{{\rm i}}}+~{{Q}_{{\rm e},{\rm i}}} \right){{V}^{\prime }}}{\rm d}x}$ and ${{f}_{\tau ,{\rm rot}}}$ is an adjustable factor of order of one.

Usually, the edge value of plasma rotation is far from being null. We use a simple model to determine it assuming there is no friction with the SoL and the only damping term is the friction with neutrals. Additionally, we assume that momentum is purely convective at the edge. The edge rotation reads:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{R}_{{\rm tot},{\rm edge}}}=\frac{\left({{S}_{\varphi ,{\rm NBI}}}+{{S}_{\varphi ,{\rm self}}}+~{{S}_{\varphi ,{\rm RF}}}+{{S}_{\varphi ,{{E}_{\parallel }}}}+~{{S}_{\varphi ,{\rm ripple}}} \right)}{\frac{\mathop{\sum }_{k\in \left\{{\rm species} \right\}}{{A}_{k}}{{n}_{k}}}{\mathop{\sum }_{k\in \left\{{\rm species} \right\}}{{Z}_{k}}{{n}_{k}}}~{{m}_{{\rm p}}}~{{n}_{{\rm out}}}~{{R}_{{\rm ref}}}-~{{F}_{n,0}}}\nonumber \end{align} \tag{ 19 }$$

where ${{n}_{{\rm out}}}$ is the flux of ions exchange through the LCFS that contains contributions from interchange, cold neutral fuelling, pellets fuelling and neutral beam injection:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{n}_{{\rm out }~ }}=\frac{2}{5}\frac{\int_{0}^{x}{\left({{Q}_{{\rm i}}}+~{{Q}_{{\rm ei}}} \right){{V}^{\prime }}}{\rm d}x}{e~{{T}_{{\rm i},a}}}+\frac{{{S}_{0}}}{{{{\bar{Z}}}_{0}}}+{{S}_{{\rm pellet}}}~+\frac{{{P}_{{\rm NBI}}}}{e{{E}_{{\rm b}0}}}.\nonumber \end{align} \tag{ 20 }$$

The first term in RHS of equation (20) is, in most of cases, the main contribution to ${{n}_{{\rm out}}}$; others terms being negligible. Nevertheless, with only this first term, the prediction of edge rotation can become unphysical during early ramp-up phase and during some transients. To prevent this to happen, terms taking into account particle sources have been added.

The third part of this calculation consists in an estimation of the radial electric field ${{E}_{{\rm r}}}\left(x \right)$, poloidal (${{V}_{\theta ,{\rm imp}}}$) and toroidal (${{V}_{\varphi ,{\rm imp}}}$) rotation of the main impurity, in the plasma equatorial plane (Z  =  Z_ref) at the low field side. We assume that the toroidal angular rotation ${{\omega }_{\varphi }}\left(x \right)$ is homothetic to a chosen kinetic profile (usually the ion temperature, which has been found to be a good proxy for ${{V}_{\varphi ,{\rm imp}}}$ in JET NBI-dominated experiments) and preserves ${{R}_{{\rm tot}}}$ with boundary condition edge value given above. By using formula (8.22 and 13.1) from [40], weighting it with mass and density and summing on all species and flux averaging it (see [2]), we are able to compute the radial electric field from ${{\omega }_{\varphi }}\left(x \right)$, pressure gradients and poloidal rotation estimation (${{V}_{\theta ,k}}={{u}_{k}}{{B}_{\theta }}$):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{E}_{{\rm r}}}=&\,\bigg [ {{\omega }_{\varphi }}\frac{{\rm d}\Psi }{{\rm d}\rho }+\frac{1}{\mathop{\sum }_{k\in \left\{{\rm species} \right\}}{{m}_{k}}{{n}_{k}}}\sum\limits_{k\in \left\{{\rm species} \right\}}\nonumber \\ &\times\left.{\left(\frac{{{m}_{k}}}{{{Z}_{k}}}\frac{{\rm d}\left({{T}_{{\rm i}}}{{n}_{k}} \right)}{{\rm d}\rho }-{{m}_{k~}}{{n}_{k}}{{u}_{k}}F\left\langle \frac{1}{{{R}^{2}}} \right\rangle \frac{{\rm d}\psi }{{\rm d}\rho } \right)}~\bigg. \right]~\left\langle \left| \nabla \rho \right| \right\rangle .\nonumber \end{align} \tag{ 21 }$$

The poloidal rotation (${{V}_{\theta ,k}}$) is estimated using either:

• The Kim formulation [37] that provides the poloidal rotation in a plasma composed of one main ion and one impurity
• Assuming no friction between ion species and using the formulation from [41]
• Assuming all ion species have the same toroidal rotation
• Assuming solid plasma rotation, i.e. ${{V}_{\theta ,k}}=0$

To compute the toroidal rotation (${{V}_{\varphi ,{\rm imp}}}$), we applied equation (14) in the [2].

METIS includes fast solvers for the computation of the source terms for the current diffusion equation (1) and the heat transport equation (4a), (4b), where electron and ion energy sources are distinguished. A source term is a radial profile which is also time-dependent and evaluated at each time slice of the METIS simulation. The source terms are also recalculated at each step of the global convergence scheme in order to obtain a self-consistent solution. Therefore, the source modules are chosen to be fast and simple enough to avoid impeding the overall computation speed of METIS. Particle and toroidal momentum profiles being not calculated from the usual transport equations, they do not require the calculation of 'source terms' in the classical definition: the external input of particles and toroidal momentum are treated differently (see sections 2.4, 2.6 and 2.7).

The source terms of the transport equations can be either prescribed as time-dependent radial profiles or calculated consistently with the simulated plasma characteristics. The source modules calculate also (when relevant) the neutron rate due to fusion reactions (either from thermal or energetic ions), as well as the fast particle pressure (parallel and perpendicular to the magnetic field) due to the heating scheme. Those quantities are not directly involved in the transport equations, though the fast particle pressure contributes to the total pressure and thus enters the equilibrium calculation as an input. They are however very useful quantities for comparison of the simulation to an experiment and data consistency applications (neutron production, stored energy) as shown for JET [42].

Several source terms have been included in METIS such as NBI, ICRH, lower hybrid (LH), electron cyclotron resonance heating (ECRF), pellets, fusion reactions of thermal ions, ripple effects, ionisation from wall recycling/gas puffing, radiative processes. In general, a simplified description is used for each source term which, on the other hand, is sufficient to have a first evaluation on their impact on the main plasma characteristics. A detailed description of the different sources is given in the appendix.

A key challenge in plasma physics is the strong non-linear coupling between plasma profiles, transport coefficients and source terms of transport equations, as sketched on figure 2. In classical 1.5D integrated modelling codes, in order to take accurately into account these non-linear couplings, a convergence loop is usually carried out at each time step, which is very expensive in terms of CPU time. In METIS, in order to use large time steps and aiming at a fast Plasma Simulator, an original numerical algorithm has been developed, similar to a waveform relaxation algorithm [43]. It consists on a global convergence loop on the whole time evolution of all plasma quantities, noted g_N at iteration N. At each iteration of this main loop, a new time evolution of all plasma quantities F(g_N) is calculated in a step-by-step update (starting from the g_N time evolution) of the quantities involved in the transport equations. This new time evolution F(g_N) is then combined to the one of the previous iteration via an oscillation damping scheme (to ease the global convergence), providing the resulting time evolution at iteration N  +  1:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{g}_{N+1}}={{\alpha }_{N}}~F\left({{g}_{N}} \right)+\left(1-{{\alpha }_{N}} \right){{g}_{N}} \nonumber \end{align} \tag{ 22 }$$

where ${{\alpha }_{N}}$ is an oscillation damping coefficient varying with the iteration number. We use ${{\alpha }_{N+1}}=0.95{{\alpha }_{N}}$ for the 21 first iterations and ${{\alpha }_{N+1}}=0.3~{{\alpha }_{N}}$ afterwards. The maximum number of iterations is set to 31. The loop terminates when the relative change between $F\left({{g}_{N}} \right)$ and g_N falls below a requested value (10⁻² in 'coarse' computation mode and 10⁻³ in 'detailed' computation mode).

Zoom In Zoom Out Reset image size

Additionally, to be able to solve ODE and PDE for any time step size, METIS uses, for time integration, an exponential integrator solver [44, 45].

The METIS code can be used from command line or through a graphical user interface (GUI). This GUI (figure 2) makes the code user-friendly and allows to parametrize and run simulations in an easy way. The GUI provides also predefined graphics and a generic tool to browse and plot data (figure 3).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

When doing simulations based on existing experiments, several plasma measurements are used to prepare the input simulation file. The access to experimental databases is handled by tools included in METIS and usable through the METIS GUI. The database access is however decoupled from the simulation itself: it is done entirely before the run. After database access, any signal can be edited and modified using the METIS GUI.

The access to experimental databases is done via a tokamak dependent routine which maps the database variables onto the METIS data structure. For some of the tokamaks already coupled to METIS (JET, TCV) database access is done via MDS. WEST and Tore Supra data access is done using the local TSLib library and the ITER integrated modelling analysis suite (IMAS) [46]. METIS is also integrated in the integrated tokamak modelling (ITM) suite of codes [8]. METIS allows also for WEST and Tore Supra to perform pre-shot simulation by means of direct access to the plasma control system (PCS) configuration data. For COMPASS, METIS uses a local mechanism to read the database. For other tokamaks (DIII-D and EAST), data preparation is based on a dedicated pre-processing tool. Data access routines also apply specific default simulation parameters for each machine. Additionally, for future machines such as ITER, DEMO or JT-60SA, METIS GUI includes a dedicated scenario generator allowing an easy preparation of scenario template. METIS is also delivered with a set of reference simulations for JT-60SA and for WEST.

METIS input data are made of the following waveforms:

• boundary condition for the current diffusion equations (plasma current or poloidal flux at LCFS)
• injected power for heat sources
• effective charge
• line-averaged density
• plasma geometry (either moments or LCFS points).
• isotopic plasma composition ($\frac{{{n}_{{\rm T}}}}{{{n}_{{\rm d}}}}$) for reactor plasma
• confinement enhancement factor
• isotopic composition of NBI ($\frac{{{\Gamma }_{{\rm T}}}}{{{\Gamma }_{{\rm D}}}}$ for reactor plasma or $\frac{{{\Gamma }_{{\rm H}}}}{{{\Gamma }_{{\rm D}}}}$ for standard plasma)

and scalar parameters describing internal physics model configurations, sources configurations and numerical scheme configuration (about 200 parameters).

Additionally, METIS internal data can be constrained by experimental data or data computed by other codes. METIS can use external data for: electron density, electron and ion temperatures, toroidal rotation, any additional heating (the shape of deposition are renormalized to the prescribed input waveform), line radiation, Z_eff (profile) and runaway electrons current.

Originally METIS has been conceived as a module of the CRONOS [2] suite of codes. Later METIS became an autonomous code and was further developed. The present version still preserves the link with CRONOS that allows preparing a METIS simulation using CRONOS data, comparing METIS results to CRONOS results and converting a METIS data set into a complete CRONOS data set. All this can be performed through the GUI of METIS and CRONOS. Additionally, METIS can use CRONOS data (for kinetic profiles and sources) instead of his internally computed data.

METIS is also completed link with the IMAS infrastructure (and on the same way with WPCD EUROfusion infrastructure, that is very similar to IMAS one). METIS can read in input simulation through UAL (Universal Access Layer) and wrote is results in IDSs (equilibrium data, kinetic profiles, all sources and many other data). METIS can be run as a standalone code in IMAS infrastructure or as an embedded actor inside a KEPLER actor.

A detailed description of the main source modules in METIS is provided in the following sections.

The plasma neoclassical resistivity and bootstrap current profiles are computed using the Sauter model [10]. Optionally, the total bootstrap current can be normalized to the Hoang scaling law [71]—the shape of the bootstrap current density profile remaining as given by the Sauter model. These quantities are used in the current diffusion equation (1).

The neutral beam injection (NBI) is described in METIS by a beam attenuation equation applied in a simplified geometry in order to calculate the fast ion source, then using an analytical solution of the Fokker–Planck equation for the fast ion distribution function.

The beam attenuation is computed for a few sub-beams spread around the tangential radius and the vertical tilt value to take into account the beam geometry. The vertical tilt is accounted for by projecting field values (i.e. density or temperature) on a tilted line. The beam intensity ($\Upsilon$) damping equation along the beam path is:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{{\rm d}\Upsilon (l)}{{\rm d}l}=-{{n}_{{\rm e}}}\left(l \right){{\sigma }_{{\rm eff}}}\left(l \right)\Upsilon \left(l \right)~\nonumber \end{align} \tag{ A.1 }$$

where $l$ is the coordinate along the beam path with the initial condition:

$\Upsilon \left(l=0 \right)={{\Upsilon }_{{\rm O}}}$ at the entrance of the plasma

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\sigma }_{{\rm eff}}}=&\,{{\sigma }_{0}}\left({{A}_{{\rm b}}},{{E}_{{\rm b}0}},{{n}_{{\rm e}}}\left(l \right),{{T}_{{\rm e}}}\left(l \right),{{Z}_{{\rm eff}}}\left(l \right) \right)\nonumber \\ &+{{\sigma }_{{\rm bp}}}\left({{A}_{{\rm b}}},{{E}_{{\rm b}0}},{{n}_{{\rm e}}}\left(l \right),{{T}_{{\rm e}}}\left(l \right),{{\xi }_{{\rm b}}}\left(l \right),{{S}_{{\rm NBI}}}\left(l \right) \right)\nonumber \end{align} \tag{ A.2 }$$

where ${{\sigma }_{0}}$ is the stopping cross section [72] and ${{\sigma }_{{\rm bp}}}$ is the increment of the stopping cross section due to fast ions [73]. The stopping cross section along the neutral path depends on the beam ion mass (${{A}_{{\rm b}}}$), the initial beam energy (${{E}_{{\rm b}0}}$), the pitch angle at the point where the neutral particles are ionized (${{\xi }_{{\rm b}}}$), the fast ion source (${{S}_{{\rm NBI}}}$), the electron density ${{n}_{{\rm e}}}$, the electron temperature ${{T}_{{\rm e}}}$ and the plasma effective charge ${{Z}_{{\rm eff}}}$.

The neutral beam path entrance point is taken in the equatorial mid-plane of the plasma (at $Z={{Z}_{0}}$) on the low field side and the radius of tangency is prescribed. The final value of Υ gives the fraction of the power that is not deposited in the plasma (shine-through).

From the power deposition (${{p}_{{\rm b}}}\left(x \right)$), we subtract the first orbit losses that are computed with a simplified model: we suppose that most of the fast ions are trapped near the plasma edge and we compute the orbit width (${{\delta }_{{\rm o}}}\left(x \right)$) as done in [74]. The fast ions created at a distance smaller than ${{\delta }_{{\rm o}}}\left(x \right)$ from the edge are lost. In order to simulate the broadening of the deposition profile due to orbit width we convolve the profile by a step function of width ${{\delta }_{{\rm o}}}\left(x \right)$.

From the power deposition (${{p}_{{\rm b}}}\left(x \right)$), we compute the fraction of the power that heats the main plasma ions using the formulae 5.4.12 in [75]. The current source associated with NBI is computed by an analytical solution of the Fokker–Planck equation in which both trapping effects and energy diffusion are neglected:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle J_{{\rm NBI}}^{{\rm fast}}=e\frac{{{p}_{{\rm b}}}}{{{E}_{{\rm b}}}}{{\tau }_{{\rm s}}}{{\xi }_{{\rm b}}}{{\left(\frac{v_{0}^{3}+v_{{\rm c}}^{3}}{v_{0}^{3}} \right)}^{\frac{2v_{\gamma }^{3}}{3v_{{\rm c}}^{3}}}}\int_{0}^{\frac{{{v}_{0}}}{{{v}_{{\rm c}}}}}{{{\left(\frac{{{z}^{3}}}{{{z}^{3}}+1} \right)}^{\frac{2v_{\gamma }^{3}}{3v_{{\rm c}}^{3}}+1}}}{\rm d}z\nonumber \end{align} \tag{ A.3 }$$

where ${{\tau }_{{\rm s}}}$ is the slowing down time, ${{v}_{{\rm c}}}$ and ${{v}_{\gamma }}$ are critical velocities and ${{v}_{0}}$ the fast ions initial velocity.

The electron back-current is computed either with the formulation found in [76] or an alternative formulation that takes into account the collisionality effect [77]. This back current is subtracted from the total NBI current to obtain the current drive by beam injection.

Ion cyclotron (IC) waves can be used in a variety of heating schemes. The user must prescribe in METIS which heating scheme is used (this scheme is applied to the whole simulation). We describe below the models used for each IC heating scheme.

A.4.1. ICRF minority heating

A.4.2. ICRH in fast wave electron heating and current drive scheme

In the fast wave (FW) electron heating and current drive scheme, it is assumed that the ICRH power heats only the electrons. The current drive efficiency is determined by a fit of experimental data [80]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{I}_{{\rm fwcd}}}={{s}_{{\rm fwcd}}}\frac{{{P}_{{\rm icrh}}}\left(t \right){{\eta }_{{\rm fwcd}}}\left(t \right)}{\frac{{{R}_{{\rm ref}}}\left(t \right)\left({{n}_{{\rm e}}}\left(t,x=0 \right) \right)}{{{10}^{20}}}}\nonumber \end{align} \tag{ A.4 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\eta }_{{\rm fwcd}}}\left(t \right)=\frac{\left(0.008\frac{{{T}_{{\rm e}}}\left(t,x=0 \right)}{{{10}^{3}}}+0.0021 \right)6}{5+{{Z}_{{\rm eff}}}\left(t \right)}\nonumber \end{align} \tag{ A.5 }$$

and where ${{s}_{{\rm fwcd}}}$ is related to the sign of the current:

${{s}_{{\rm fwcd}}}=1$ for co-current (i.e. in the direction of toroidal plasma current), ${{s}_{{\rm fwcd}}}=-1$ for counter current and ${{s}_{{\rm fwcd}}}=0$ for FW electron heating scheme.

The shape of the power deposition profile is:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{p}_{{\rm fw}}}\propto \frac{{{n}_{{\rm e}}}\left(t,x \right){{T}_{{\rm e}}}\left(t,x \right)}{{{B}_{{\rm out}}}{{\left(t,x \right)}^{3}}~{{V}^{\prime }}\left(t,x \right)} \nonumber \end{align} \tag{ A.6 }$$

where ${{B}_{{\rm out}}}\left(t,x \right)~$ is the total magnetic field in the middle plane at the low magnetic field side.

with ${{\rho }_{{\rm m}}}\int_{0}^{1}{{{p}_{{\rm fw}}}\left(t,x \right){{V}^{\prime }}\left(t,x \right){\rm d}x={{P}_{{\rm ICRH}}}\left(t \right)}$

The shape is for the current density profile is:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{j}_{{\rm fwcd}}}\propto \frac{{{n}_{{\rm e}}}\left(t,x \right){{T}_{{\rm e}}}\left(t,x \right){{\sigma }_{\parallel }}\left(t,x \right)}{{{B}_{{\rm out}}}{{\left(t,x \right)}^{3}}~{{V}^{\prime }}\left(t,x \right)}~\nonumber \end{align} \tag{ A.7 }$$

that is normalised to the total current ${{I}_{{\rm fwcd}}}$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\rho }_{{\rm m}}}\int_{0}^{1}{{{j}_{{\rm fwcd}}}\left(t,x \right){{S}^{\prime }}\left(t,x \right){\rm d}x={{I}_{{\rm fwcd}}}\left(t \right)}\nonumber \end{align}$$

The determination of the power deposition of EC waves requires a ray tracing code, but this type of code would be too slow to be used in METIS. Therefore, the location of the maximum of the power deposition profile, x_eccd, has to be prescribed by the user. The shape of the power deposition is a Gaussian curve with a width, δ_eccd, determined from the width of the resonance [75], with constants adjusted to fit the results of the REMA ray-tracing code [81]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{p}_{{\rm eccd}}}\left(x,t \right)={{p}_{{\rm eccd,0}}}{{e}^{\frac{{{\left(x-{{x}_{{\rm eccd}}}\left(t \right) \right)}^{2}}}{2{{\delta }_{{\rm eccd}}}\left(t \right)}}}\nonumber \end{align} \tag{ A.8 }$$

where ${{\delta }_{{\rm eccd}}}\left(t \right)=\sqrt{\frac{1}{4}{{\left(\frac{{{v}_{{\rm th}}}}{c} \right)}^{4}}+{{\alpha }_{{\rm eccd}}}\frac{{{I}_{{\rm eccd}}}}{{{P}_{{\rm eccd}}}}{{\left(\frac{{{v}_{{\rm th}}}}{c} \right)}^{2}}}$ with ${{v}_{{\rm th}}}=\sqrt{\frac{2e{{T}_{{\rm e}}}\left(t,{{x}_{{\rm eccd}}} \right)}{{{m}_{{\rm e}}}}}$, ${{\alpha }_{{\rm eccd}}}\simeq 1$ and ${{p}_{{\rm eccd},0}}$ satisfying ${{\rho }_{{\rm m}}}\int_{0}^{1}{{{p}_{{\rm eccd}}}\left(t,x \right){{V}^{\prime }}\left(t,x \right)}$${{\rm d}x={{P}_{{\rm eccd}}}\left(t \right)}$.

The current source profile ${{j}_{{\rm eccd}}}\left(t,x \right)$ has the same shape as the power deposition profile. The total current is computed using a simple scaling [81]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle& {{\rho }_{{\rm m}}}\int_{0}^{1}{{{j}_{{\rm eccd}}}\left(t,x \right){{S}^{\prime }}\left(t,x \right){\rm d}x={{I}_{{\rm eccd}}}\left(t \right)\,{\rm with}}\nonumber \\ &{{{I}_{{\rm eccd}}}}{\left(t \right)={{s}_{{\rm eccd}}}\frac{{{\Gamma }_{{\rm LH}\in {\rm ECCD}}}{{\eta }_{{\rm eccd}}}\left(t \right){{P}_{{\rm eccd}}}\left(t \right)}{{{n}_{{\rm e}}}\left(t,{{x}_{{\rm eccd}}}\left(t \right) \right){{R}_{{\rm ref}}}\left(t \right)}}\nonumber \end{align} \tag{ A.9 }$$

where ${{s}_{{\rm eccd}}}$ is the direction of the wave injection:

${{s}_{{\rm eccd}}}=1$ for co-current, ${{s}_{{\rm eccd}}}=-1$ for counter current and ${{s}_{{\rm eccd}}}=0$ for normal injection.

and where ${{\Gamma }_{{\rm LH}\in {\rm ECCD}}}$ is the synergy factor with lower hybrid (equal to 1 by default) and where:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\eta }_{{\rm eccd}}}=\frac{\frac{1}{1+\frac{100}{\frac{{{T}_{{\rm e}}}\left(t,{{x}_{{\rm eccd}}} \right)}{{{10}^{3}}}}\left[ 1-\left(1+\frac{5+{{Z}_{{\rm eff}}}}{3\left(1+{{Z}_{{\rm eff}}} \right)} \right){{\left(\sqrt{2}{{\mu }_{{\rm t}}} \right)}^{\frac{5+{{Z}_{{\rm eff}}}}{1+{{Z}_{{\rm eff}}}}}} \right]6}}{5+{{Z}_{{\rm ef\!f}}}}\nonumber \end{align} \tag{ A.10 }$$

with ${{\mu }_{{\rm t}}}=\sqrt{\frac{{{a}_{{\rm ref}}}{{x}_{{\rm eccd}}}\left(1+\cos \left({{\theta }_{{\rm pol}}} \right) \right)}{{{R}_{{\rm ref}}}+{{a}_{{\rm ref}}}{{x}_{{\rm eccd}}}\cos \left({{\theta }_{{\rm pol}}} \right)}}$.

In some configurations, a synergy between ECCD and LHCD may increase the effect of the EC power on the total driven current. The factor ${{\Gamma }_{{\rm LH}\in {\rm ECCD}}}$ accounts for this effect. This factor must be prescribed except for Tore Supra [82], for which the following expression can be applied:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle& {{s}_{{\rm eccd}}}{{\Gamma }_{{\rm LH}\in {\rm ECCD}}}=~{{s}_{{\rm eccd}}}+{{H}_{{\rm MHD}}}\frac{\min \left(\frac{2}{3}{{P}_{{\rm LH}}},{{P}_{{\rm ECCD}}} \right)}{\max \left(1,{{P}_{{\rm ECCD}}} \right)}\nonumber \\ &\frac{{{s}_{{\rm ECCD}}}+1}{2}\max \left(0,\frac{{{\eta }_{{\rm LH}}}}{{{\eta }_{{\rm eccd}}}} \right)\sqrt{\frac{{{\left(1-{{x}_{{\rm eccd}}}^{2} \right)}^{{{\nu }_{{\rm T}}}}}}{2}}\nonumber \end{align} \tag{ A.11 }$$

We designed a heuristic model for Lower Hybrid power deposition and current drive that is based on observation in present devices and on prediction made by the code C3PO/LUKE [66]. The model for the power deposition is purely phenomenological. The current drive efficiency is based on a simple physical assumption with one free parameter tuned with the help of experimental measurements and code predictions. In METIS, a positive parallel refractive index n_|| generates a co-current source and a negative n_|| generates a counter-current source, whatever is the real geometry and orientation of the magnetic field and plasma current in the tokamak. The computation is done for the main positive n_|| lobe of the spectra that generates a co-current source and for the main negative n_|| lobe that generates a counter-current source. The power deposition model is based on a probabilistic ad-hoc formulation with penalization for the wave absorption. This model only takes into account the limits of the propagation domain and the Landau absorption criterion. The probability formulation has been obtained by a trial and error method, guided by physical knowledge of the problem. The probability function of absorption, on the flux surface labelled per x (varying from 0 on magnetic axis to 1 on LCFS), reads:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{P}_{{\rm abs}}}\left(x \right)=~C{{P}_{{\rm Landau }~ }}\left(x \right)~{{L}_{{\rm acc}}}\left(x \right)~{{L}_{{\rm lc}}}\left(x \right)~{{L}_{{\rm hc}}}\left(x \right)\nonumber \end{align} \tag{ A.12 }$$

with C is normalization constant.

The probability of absorption by Landau effect reads:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{P}_{{\rm Landau}}}\left(x \right)=~{{{\rm e}}^{-{{\left(\frac{{{n}_{{\rm Landau}}}\left(x \right)-~{{n}_{\parallel ,0}}-\frac{\Delta {{n}_{\parallel }}\left(x \right)}{2}}{\Delta {{n}_{\parallel }}\left(x \right)} \right)}^{2}}}}.\nonumber \end{align} \tag{ A.13 }$$

The penalization of absorption probability due to density accessibility limit reads:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{L}_{{\rm acc}}}\left(x \right)=\min \left(1,{{{\rm e}}^{\frac{{{n}_{{\rm Landau}}}\left(x \right)-{{n}_{||,{\rm acc}}}\left(x \right)-\Delta {{n}_{\parallel ,0}}}{\Delta {{n}_{\parallel ,0}}}}} \right).\nonumber \end{align} \tag{ A.14 }$$

The penalization of absorption probability due to lower n_|| possible propagation reads:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{L}_{{\rm lc}}}\left(x \right)=\min \left(1,{{{\rm e}}^{\frac{{{n}_{{\rm Landau}}}\left(x \right)-~\left({{n}_{\parallel ,{\rm lc}}}\left(x \right)+\frac{\Delta {{n}_{\parallel }}\left(x \right)}{2} \right)}{\Delta {{n}_{\parallel }}\left(x \right)}}} \right).\nonumber \end{align} \tag{ A.15 }$$

The penalization of absorption probability due to higher n_|| possible propagation reads:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{L}_{{\rm hc}}}\left(x \right)=\min \left(1,{{{\rm e}}^{-\frac{{{n}_{{\rm Landau}}}\left(x \right)-~\left({{n}_{\parallel ,{\rm hc}}}\left(x \right)-\frac{\Delta {{n}_{\parallel }}\left(x \right)}{2} \right)}{\Delta {{n}_{\parallel }}\left(x \right)}}} \right)\nonumber \end{align} \tag{ A.16 }$$

where the ad-hoc broadening of the wave spectrum is:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Delta {{n}_{\parallel }}\left(x \right)=\left(\Delta {{n}_{\parallel ,0}}+{{n}_{\parallel ,0}}\frac{{{\Delta }_{{\rm a}}}}{{{R}_{{\rm axe}}}\left(x \right)} \right)\left(1+\frac{{{n}_{{\rm Landau}}}\left(x \right)}{{{n}_{\parallel ,0}}} \right).\nonumber \end{align} \tag{ A.17 }$$

In which $\Delta {{n}_{\parallel ,0}}=\frac{c}{{{f}_{{\rm LH}}}{{\Delta }_{{\rm a}}}~}$ is the spectrum width of the launched wave.

The Landau resonance is taken as the simple following expression:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{n}_{{\rm Landau}}}\left(x \right)=~\frac{6.5}{\sqrt{{{10}^{-3}}~{{T}_{{\rm e}}}\left(x \right)}}.\nonumber \end{align} \tag{ A.18 }$$

The accessibility limit reads:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{n}_{\parallel ,{\rm acc}}}\left(x \right)=~\frac{{{\omega }_{{\rm pe}}}\left(x \right)}{{{\omega }_{{\rm ce}}}\left(x \right)}+~\sqrt{1+\frac{\omega _{{\rm pe}}^{2}\left(x \right)}{\omega _{{\rm ce}}^{2}\left(x \right)}-\frac{\omega _{{\rm pi}}^{2}\left(x \right)}{{{\omega }^{2}}}}.\nonumber \end{align} \tag{ A.19 }$$

The lower bound of propagation domain reads:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{n}_{\parallel ,{\rm lc}}}\left(x \right)=~\frac{{{n}_{\parallel ,0}}}{1+\frac{\rho \left(x \right)}{q\left(x \right){{R}_{{\rm axe}}}\left(x \right)~}\sqrt{-\frac{p\left(x \right)}{s\left(x \right)}}}\nonumber \end{align} \tag{ A.20 }$$

The upper bound of propagation domain reads:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{n}_{\parallel ,{\rm hc}}}\left(x \right)=~\frac{{{n}_{\parallel ,0}}}{1-\frac{\rho \left(x \right)}{q\left(x \right){{R}_{{\rm axe}}}\left(x \right)~}\sqrt{-\frac{p\left(x \right)}{s\left(x \right)}}} \nonumber \end{align} \tag{ A.21 }$$

where $\rho \left(x \right)=~\sqrt{\frac{\Phi \left(x \right)}{{{\Phi }_{\max }}}}$ and $\Phi \left(x \right)$ is the toroidal flux (Wb), $q\left(x \right)$ the safety factor, ${{R}_{{\rm axe}}}\left(x \right)$ the geometrical axes of each flux surface (m), ${{n}_{\parallel ,0}}$ is the absolute value of the launched wave parallel index of refraction, $\omega =2\pi ~{{f}_{{\rm LH}}}$, is the frequency of the LH wave (${{H}_{z}}$), ${{\Delta }_{{\rm a}}}$ is the toroidal width of the active part of the launcher (m) and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle p\left(x \right)=1-\frac{\omega _{{\rm pe}}^{2}\left(x \right)}{{{\omega }^{2}}}-\frac{\omega _{{\rm pi}}^{2}\left(x \right)}{{{\omega }^{2}}}\nonumber \end{align} \tag{ A.22 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle s\left(x \right)=1+\frac{\omega _{{\rm pe}}^{2}\left(x \right)}{\omega _{{\rm ce}}^{2}\left(x \right)}-\frac{\omega _{{\rm pi}}^{2}\left(x \right)}{{{\omega }^{2}}}\nonumber \end{align} \tag{ A.23 }$$

${\rm where}\,{{\omega }_{{\rm pe}}}=\sqrt{\frac{{{e}^{2}}{{n}_{{\rm e}}}\left(x \right)}{{{m}_{{\rm e }~ }}{{\varepsilon }_{0}}}}$ is the electron plasma frequency, ${{\omega }_{{\rm pi}}}=\sqrt{\frac{{{e}^{2}}{{n}_{{\rm i}}}\left(x \right)}{A~{{m}_{{\rm p }~ }}{{\varepsilon }_{0}}}}$ is the ion plasma frequency, ${{\omega }_{{\rm ce}}}=~\frac{eB\left(x \right)}{{{m}_{{\rm e}}}}$ is electron cyclotron frequency, ${{T}_{{\rm e}}}\left(x \right)$ the electron temperature (eV), ${{n}_{{\rm e}}}\left(x \right)$ the electron density (${{{\rm m}}^{-3}}$), ${{n}_{{\rm i}}}\left(x \right)$ the main ion density (${{{\rm m}}^{-3}}$), $B(x)$ the total magnetic field averaged on a flux surface (T), e the charge of the electron (C), ${{m}_{{\rm e}}}$ the electron mass (kg), ${{\varepsilon }_{0}}$ the vacuum permittivity (${\rm F}\,{{{\rm m}}^{-1}}$), A the number of mass of the main ion and c the speed of light in vacuum (${\rm m}\,{{{\rm s}}^{-1}}$). The constant C is computed to have ${{\rho }_{\max }}\int_{0}^{1}{P\left(x \right){{V}^{\prime }}\left(x \right){\rm d}x=1}$. If the input power is ${{P}_{{\rm in}}}$, the local power deposition is ${{P}_{{\rm tot}}}\left(x \right)={{P}_{{\rm in}}}~P\left(x \right)$.

We now describe how the amount of current drive by LH wave is computed. We restart from the Fisch formulation of the current drive efficiency [83] and we add a correction to take into account the quasi linear effect ${{Q}_{{\rm l}}}\left(x \right)$ and a penalization term to take into account the accessibility $\newcommand{\e}{{\rm e}} {{\eta }_{{\rm acc}}}$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\eta }_{{\rm LHCD}}}\left(x \right)=~{{\eta }_{0}}\left(x \right)~{{\eta }_{{\rm acc}}}~{{Q}_{{\rm l}}}\left(x \right)~\frac{\omega _{2}^{2}-\omega _{1}^{2}}{\ln \left(\frac{{{\omega }_{2}}}{{{\omega }_{1}}} \right)}.\nonumber \end{align} \tag{ A.24 }$$

The plasma current density reads:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{j}_{{\rm LHCD}}}\left(x \right)={{\eta }_{{\rm LHCD}}}\left(x \right)\frac{p\left(x \right)}{{{n}_{{\rm e}}}\left(x \right)}.\nonumber \end{align} \tag{ A.25 }$$

With ${{\omega }_{1}}=\frac{1}{{{n}_{{\rm Landau}}}\left(x \right)}$, ${{\omega }_{2}}=\frac{1}{{{n}_{\parallel ,0}}}$, ${{Q}_{{\rm l}}}\left(x \right)=\frac{3}{2}+~{\scriptscriptstyle 1/{}_2}\tanh \left[ \ln \left(10~{{D}_{\parallel }}\left(x \right) \right) \right]$, $\newcommand{\e}{{\rm e}} ~{{\eta }_{{\rm acc}}}=\min \left(1,{{e}^{\frac{{{n}_{\parallel ,0}}-{{n}_{\parallel ,{\rm acc}}}\left(1 \right)~}{\Delta {{n}_{\parallel ,0}}}}} \right)$ and $\newcommand{\e}{{\rm e}} {{\eta }_{0}}\left(x \right)=\frac{{{\eta }_{{\rm m}}}}{\ln {{\Lambda }_{{\rm ee}}}\left(x \right)}$ [84], where $\ln {{\Lambda }_{{\rm ee}}}\left(x \right)\sim 14.9+\frac{1}{2}\ln \left(\frac{{{n}_{{\rm e}}}\left(x \right)}{{{10}^{20}}} \right)+\ln \left(\frac{{{T}_{{\rm e}}}\left(x \right)}{{{10}^{3}}} \right)$ is the Coulomb logarithm and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{D}_{\parallel }}\left(x \right)=0.32~\frac{{{P}_{{\rm in}}}}{{{10}^{6}}}~\frac{n_{{\rm Landau}}^{2}\left(x \right)\sqrt{\frac{{{T}_{{\rm e}}}\left(x \right)}{{{10}^{3}}}}}{\sqrt{s\left(x \right)}\Delta {{n}_{\parallel }}\left(x \right)\max \left({{\Delta }_{{\rm a}}},\rho \left(x \right) \right)~{{R}_{{\rm axe}}}\left(x \right){{\left(\frac{{{n}_{{\rm e}}}\left(x \right)}{{{10}^{19}}} \right)}^{\frac{3}{2}}}}\nonumber \end{align} \tag{ A.26 }$$

is a proxy of the quasi linear diffusion coefficient.

The constant $\newcommand{\e}{{\rm e}} {{\eta }_{{\rm m}}}$ is set to $3.1~{{10}^{21}}$ to fit experimental results and simulations.

If ${{Z}_{{\rm eff }~ }}>1$, the effect on current drive efficiency must be taken into account. The effect of trapping particles has been included by an analytical formula [85]. Finally the net current drive is:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{j}_{{\rm LHCD}}}\left(x \right)={{\eta }_{{\rm LHCD}}}\left(x \right)\frac{p\left(x \right)}{{{n}_{{\rm e}}}\left(x \right)}~\frac{1-~{{\varepsilon }^{\frac{5+~{{Z}_{{\rm eff}}}\left(x \right)}{2\left(1+~{{Z}_{{\rm eff}}}\left(x \right) \right)}}}}{5+{{Z}_{{\rm eff}}}\left(x \right)}~\nonumber \end{align} \tag{ A.27 }$$

with $\varepsilon =\frac{a\left(x \right)}{{{R}_{{\rm axe}}}\left(x \right)}$

The fusion reactions between deuterium and tritium, between deuterium and deuterium and tritium and tritium are calculated, including the thermal plasma, beam-plasma and beam-beam reactions. The tritium-tritium reaction is negligible compared to the rate of D–T reactions in a fusion reactor but it is important for the neutron diagnostic.

For each of the considered fusion reactions, METIS calculates the resulting neutron production rate, power deposition profiles on the plasma electrons and ions, as well as (for D–T reactions) the average Helium ash production rate, which is used to calculate the Helium average density following the procedure described in section D.1. The fast alpha particles originating from the D–T reactions also drive a bootstrap-like current, which is estimated by METIS.

Reactivities, i.e. $\left\langle \sigma v \right\rangle$ and cross sections used in METIS are those of [86] but for tritium–tritium cross section library [87] is used.

A.8.1. Thermal reactions

The fusion power source is simply computed using deuterium and tritium density profiles and the ion temperature profile:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{S}_{\alpha ,{\rm th}}}={{n}_{{\rm D}}}\left(t,x \right){{n}_{{\rm T}}}\left(t,x \right){{\left\langle \sigma v \right\rangle }_{{\rm T}\left({\rm d},{\rm n} \right){\rm He}4}}\left({{T}_{i}}\left(t,x \right) \right)\nonumber \end{align} \tag{ A.28 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{p}_{\alpha ,{\rm th}}}=e{{E}_{\alpha ,{\rm T}\left({\rm d},{\rm n} \right){\rm He}4}}{{S}_{\alpha ,{\rm th}}}\nonumber \end{align} \tag{ A.29 }$$

where ${{E}_{\alpha ,{\rm T}\left(d,n \right){\rm He}4}}=3.56\times {{10}^{6}}$ eV

For DD reactions there are two channels:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{S}_{{\rm p},{\rm DD},{\rm th}}}=\frac{1}{2}{{n}_{{\rm D}}}{{\left(t,x \right)}^{2}}{{\left\langle \sigma v \right\rangle }_{{\rm D}\left({\rm d},{\rm p} \right)T}}\left({{T}_{{\rm i}}}\left(t,x \right) \right)\nonumber \end{align} \tag{ A.30 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{S}_{{\rm n},{\rm DD},{\rm th}}}=\frac{1}{2}{{n}_{{\rm D}}}{{\left(t,x \right)}^{2}}{{\left\langle \sigma v \right\rangle }_{{\rm D}\left({\rm d},n \right){\rm He}3}}\left({{T}_{{\rm i}}}\left(t,x \right) \right).\nonumber \end{align} \tag{ A.31 }$$

The associated power source in the plasma is:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{p}_{{\rm DD},{\rm th}}}=e\left({{E}_{{\rm p},{\rm D}\left({\rm d},{\rm p} \right){\rm T}}}+{{E}_{{\rm T},{\rm D}\left({\rm d},{\rm p} \right){\rm T}}} \right){{S}_{{\rm p},{\rm DD},{\rm th}}}+e{{E}_{{\rm He}3,{\rm D}\left({\rm d},{\rm n} \right){\rm He}3}}{{S}_{{\rm n},{\rm DD},{\rm th}}}\nonumber \end{align} \tag{ A.33 }$$

where ${{E}_{{\rm p},{\rm D}\left({\rm d},{\rm p} \right){\rm T}}}=3.02\times {{10}^{6}}{\rm eV}$; ${{E}_{{\rm T},{\rm D}\left({\rm d},{\rm p} \right){\rm T}}}=1.01\times {{10}^{6}}{\rm eV}$ and ${{E}_{{\rm He}3,{\rm D}\left({\rm d},{\rm n} \right){\rm He}3}}=0.82\times {{10}^{6}}\,{\rm eV}$

From the power deposition, as in the NBI case, we compute the fraction of the power that heats the main plasma ions using the formulae 5.4.12 in [75].

A.8.2. Neutron production rate

A.9.1. Introduction

A.9.2. Line radiation and Bremsstrahlung

The radiative power is computed from the temperature and density profiles of each species using the radiative collisional equilibrium [89]. This gives an estimate of the profile of power radiated by line transition and thermal Bremsstrahlung by using the cooling rate l(Z_k, T_e) that combines line and Bremsstrahlung radiation (cooling rates are extracted from the ADAS data base [90]):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{p}_{{\rm rad}}}\left(t,x \right)={{n}_{{\rm e}}}\left(t,x \right)\underset{k\in \left\{{\rm species} \right\}}{\mathop \sum }\,l\left({{Z}_{k}},{{T}_{{\rm e}}}\left(t,x \right) \right){{n}_{k}}\left(t,x \right).\nonumber \end{align} \tag{ A.34 }$$

This radiative power profile is integrated over the plasma volume to compute the total radiated power ${{P}_{{\rm rad}}}={{\rho }_{{\rm m}}}\int_{0}^{1}{{{V}^{\prime }}{{p}_{{\rm rad}}}}{\rm d}x$.

The total radiated power due to line radiation and bremsstrahlung can optionally be normalized to the Matthews law [91, 92]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{P}_{{\rm rad},{\rm mat}}}=\frac{{{10}^{6}}}{4.5~\left({{Z}_{{\rm eff}}}-{{Z}_{{\rm main}}} \right){{\left(\frac{{\bar{n}}}{{{10}^{20}}} \right)}^{1.89}}}{{S}^{0.94}}Z_{\max }^{-0.12} \nonumber \end{align} \tag{ A.35 }$$

where $S$ is the external surface of the confined plasma (i.e. the surface of the LCFS) and ${{Z}_{{\rm main}}}$ is the charge number of main plasma ions (${{Z}_{{\rm main}}}=1$ for hydrogen, deuterium or deuterium/tritium plasma and ${{Z}_{{\rm main}}}=2$ for helium plasma). As there is always, at least a small amount of impurities in the plasma, we have always ${{Z}_{{\rm eff}}}>{{Z}_{{\rm main}}}$ preventing ${{P}_{{\rm rad},{\rm mat}}}$ to diverge.

To calculate the power flowing through the LCFS ${{P}_{{\rm loss}}}$, we need to separate the bremsstrahlung radiation ${{P}_{{\rm brem}}}$ (fully radiated by the confined plasma and thus fully subtracted for the ${{P}_{{\rm loss}}}$ calculation) and the line radiation ${{P}_{{\rm line}}}$ (partially radiated by the confined plasma and partially in the scrape-off layer, thus only a fraction of it is subtracted for the ${{P}_{{\rm loss}}}$ calculation). We use the following expression for the bremsstrahlung radiated power coming from volume integration of local expressions [93, 94]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle& {{p}_{{\rm brem}}}\left(t,x \right)=4.8562\times {{10}^{3}}~\frac{{{n}_{{\rm e}}}\left(t,x \right)}{{{10}^{20}}}\sqrt{\frac{{{T}_{{\rm e}}}\left(t,x \right)}{{{10}^{3}}}}\nonumber \\ &\underset{k\in \left\{{\rm species} \right\}}{\mathop \sum }\,{{C}_{{\rm Gaunt}}}\left({{T}_{{\rm e}}}\left(t,x \right),{{Z}_{k}} \right)~Z_{k}^{2}\frac{{{n}_{k}}\left(t,x \right)}{{{10}^{20}}}\nonumber \end{align} \tag{ A.36 }$$

with ${{C}_{{\rm Gaunt}}}$ the Gaunt factor that is given by a tabulated function [92] (taking ${{C}_{{\rm Gaunt}}}=1.2$ gives an accuracy within about 20%) and ${{P}_{{\rm brem}}}={{\rho }_{{\rm m}}}\int_{0}^{1}{{{V}^{\prime }}{{p}_{{\rm brem}}}}{\rm d}x$

We then define:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{P}_{{\rm line}}}={{P}_{{\rm rad}}}-{{P}_{{\rm brem}}}\nonumber \end{align} \tag{ A.37 }$$

We neglect here a small correction: low temperature finite Rydberg energy, electron–electron bremsstrahlung and re-absorption, but we introduce the high temperature corrections due to relativistic effects [95]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{P}_{{\rm brem},{\rm rel}}}={{P}_{{\rm brem}}}{{X}_{{\rm rel}}}\left(t,x \right)\nonumber \end{align} \tag{ A.38 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{X}_{{\rm rel}}}\left(t,x \right)=\left(1+\frac{2{{T}_{{\rm e}}}\left(t,x \right)}{511~{{10}^{3}}} \right)\left(1+\frac{2}{\left({{Z}_{{\rm eff}}}\left(t,x \right) \right)\left(1-\frac{1}{1+\frac{{{T}_{{\rm e}}}\left(t,x \right)}{51\,1~{{10}^{3}}}} \right)} \right).\nonumber \end{align} \tag{ A.39 }$$

With the same formulation, we estimate the total radiative loss in the SoL (P_SOL) taking an exponential decrease of the profile with a characteristic length ${{\lambda }_{{\rm SOL}}}$. This factor is either prescribed as a fraction of the minor radius or computed with a scaling law [96].

A.9.3. Cyclotron radiation

The cyclotron radiation power loss is given by the Albajar scaling [97, 98]. These expressions take into account the absorption and re-emission of the cyclotron radiation in the plasma.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle& {{P}_{{\rm cyclo}}}\left(W \right)=3.84~{{10}^{-2}}{{\left(1-{{r}_{{\rm w}}} \right)}^{1/2}}R{{a}^{1.38}}{{\kappa }^{0.79}}B_{{\rm t}}^{2.62}n_{{\rm e},{\rm c}}^{0.38}\nonumber \\ &\times {{T}_{{\rm e},{\rm c}}}{{\left(16+{{T}_{{\rm e},{\rm c}}} \right)}^{2.61}}{{\left(1+0.12\frac{{{T}_{{\rm e,c}}}}{P_{a0}^{0.14}} \right)}^{-1.51}}K\left({{\alpha }_{n}},{{\alpha }_{{\rm T}}},{{\beta }_{{\rm T}}} \right)G\left(\frac{R}{a} \right)\nonumber \end{align} \tag{ A.40 }$$

with

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{P}_{{\rm a},{\rm c}}}=6.04\times {{10}^{3}}\frac{a{{n}_{{\rm e},{\rm c}}}}{{{B}_{{\rm t}}}}\nonumber \end{align} \tag{ A.41 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle &K\left({{\alpha }_{{\rm n}}},{{\alpha }_{{\rm T}}},{{\beta }_{{\rm T}}} \right)={{\left({{\alpha }_{{\rm n}}}+3.87~{{\alpha }_{{\rm T}}}+1.46 \right)}^{-0.79~}}\nonumber \\ &\times{{\left(1.98+{{\alpha }_{{\rm T}}} \right)}^{1.38}}\beta _{{\rm T}}^{2.14}~{{\left(\beta _{{\rm T}}^{1.53}+1.87~{{\alpha }_{{\rm T}}}-0.16 \right)}^{-1.33}}\nonumber \end{align} \tag{ A.42 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle G\left(A \right)=0.93~\left(1+0.85{{e}^{-~0.82A}} \right)\nonumber \end{align} \tag{ A.43 }$$

and ${{n}_{{\rm e},{\rm c}}}=\frac{{{n}_{{\rm e}}}\left(t,x=0 \right)}{{{10}^{20}}}$, $~{{T}_{{\rm e},{\rm c}}}=\frac{{{T}_{{\rm e}}}\left(t,x=0 \right)}{{{10}^{3}}}\,{\rm and}\,{{p}_{{\rm a},{\rm c}}}={{6.0410}^{3}}{{a}_{{\rm ref}}}\frac{{{n}_{{\rm e},{\rm c}}}}{{{B}_{{\rm ref}}}}~$ where ${{\beta }_{{\rm T}}}$ and ${{\alpha }_{{\rm T}}}$ are determined from the best fit of the electron temperature profile with the shape:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{T}_{{\rm e},{\rm s}}}=\left({{T}_{{\rm e},0}}-~{{T}_{{\rm e},{\rm a}}} \right){{\left(1-{{x}^{{{\beta }_{{\rm T}}}}} \right)}^{{{\alpha }_{{\rm T}}}}}+~{{T}_{{\rm e},{\rm a}}}\,{\rm with}\,x\in [{\rm 0,0}{\rm .9 }] \nonumber \end{align} \tag{ A.44 }$$

and with ${{r}_{{\rm w}}}$ is the effective wall reflection coefficient for cyclotron radiation.

The radial profile of cyclotron radiation loss power is taken as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle& {{p}_{{\rm cyclo}}}\left(t,x \right)=~{{p}_{{\rm c}0}}~\left(t \right)~{{R}_{{\rm axis}}}\left(t,x \right)~{{K}_{x}}{{\left(t,x \right)}^{0.79}}\nonumber \\ &\times{{\left(\frac{{{F}_{{\rm dia}}}\left(t,x \right)}{~{{R}_{{\rm axis}}}\left(t,x \right)} \right)}^{2.62}}{{n}_{{\rm e}}}{{\left(t,x \right)}^{0.38}}~{{T}_{{\rm e}}}{{\left(t,x \right)}^{3.61}}\nonumber \end{align} \tag{ A.45 }$$

With ${{p}_{{\rm c}0}}$ verifying ${{p}_{{\rm c}0}}=\frac{{{P}_{{\rm cyclo}}}}{{{\rho }_{{\rm m}}}\int_{0}^{1}{{{V}^{\prime }}{{p}_{{\rm pcyclo}}}}{\rm d}x}$ where ${{K}_{x}}$ is the elongation of each flux surface and ${{F}_{{\rm dia}}}$ is the diamagnetic function, both taken from equilibrium.

We have verified that ${{p}_{{\rm cyclo}}}$ computed with this formula matches well CYTRAN [99] and EXACTEC [99] results for reactor case but not account for heat transport from center to edge as provided by EXACTEC code.

METIS can model pellet injection. The user must prescribe the maximum ${{x}_{{\rm pellet}}}$ and optionally the width ${{\Delta }_{{\rm pellet}}}$ of the pellet deposition (${{S}_{{\rm pellet}}}$), assumed to have a Gaussian shape (centered at ${{x}_{{\rm pellet}}}$ of width ${{\Delta }_{{\rm pellet}}}$), and the relative density increase provided by pellet injection ${{f}_{{\rm pellet}}}$. If ${{\Delta }_{{\rm pellet}}}$ is not provided, the shape of the deposition is assumed to follow the NGS model [100] in which the injection velocity is adjusted to obtain the desired ${{x}_{{\rm pellet}}}$ and the radius of the pellet is adjusted to have matter deposition corresponding to ${{f}_{{\rm pellet}}}\left(t \right)$.

The pellet injection does not change the line averaged density as it is supposed already included in the prescribed $\overline{{{n}_{{\rm e}}}}$ waveform (see section 2.4). Only the density profile shape is modified: a perturbation of density profile due to the pellet deposition is computed ($\Delta {{n}_{{\rm e}}}=~\left(1-x \right)~{{f}_{{\rm pellet}}}~{{S}_{{\rm pellet }~ }}$, where $\left(1-x \right)~$ accounts for efficiency of pellet deposition to change the density profile) and then is added to the original density profile (without pellets, calculated as explained in section 2.4) and the obtained density profile is then renormalized to the line averaged density keeping unchanged the edge density. The resulting density profile has a new shape with a new peaking factor.

Pellets can be described either as a continuous or discrete effect on the density profile. In the latter case, the pellet injection events are detected from sharp variations in the prescribed line averaged density waveform.

The neutral source is computed to account for fuelling with recycling particles and gas puffing. This source gives the ionisation heat source at the edge, the friction that slows down edge toroidal rotation and the electron flux across the plasma. We consider two populations. The first (cold) is made of particles coming from gas puffing and is assumed to have low temperature (typically the temperature of the vacuum vessel). The second (hot) is made of recycling particles and is assumed to have local ion temperature due to the charge exchange mechanism. Both sources are computed using a one-dimensional model based on diffusion, where the diffusion coefficient is deduced from the mean free path [101, 102]. The mean free path is computed using the charge-exchange and ionisation rate. Only monoatomic species are taken into account. Molecular dissociation and other effects are not considered. The detail of the model can be found in the [103].

The amount of hot neutrals entering the plasma is the recycling fraction (${{f}_{n0}}$) of ions flux in the SOL. This fraction depends on the configuration, i.e. if the plasma is in the limiter or the divertor configuration. The flux of ions in the SoL (sufficiently far from limiter or divertor target), is computed as ${{n}_{{\rm e},{\rm a}}}~S_{{\rm LCFS}}^{\prime }{{\delta }_{{\rm SOL}}}{{c}_{{\rm s,LCFS}}}~+~\frac{{\rm d}N}{{\rm d}t}$, given the flux of neutral entering in the LCFS:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{S}_{{\rm 0,in}}}=~{{f}_{n0~}}\left({{n}_{{\rm e},{\rm a}}}~S_{{\rm LCFS}}^{\prime }{{\delta }_{{\rm SOL}}}{{c}_{{\rm s,LCFS}}}~+~\frac{{\rm d}N}{{\rm d}t}-~{{S}_{{\rm NBI}}} \right)\left(1-~{{f}_{{\rm pellet}}} \right) \nonumber \end{align} \tag{ A.46 }$$

where ${{n}_{{\rm e},{\rm a}}}$ is the electron density at LCFS, $S_{{\rm LCFS}}^{\prime }$ is the perimeter of the poloidal plasma section, ${{\delta }_{{\rm SOL}}}$ is the SoL width, ${{c}_{{\rm s},{\rm LCFS}}}$ is the sound speed at LCFS, N is the total number of electrons in the plasma, ${{S}_{{\rm NBI}}}$ is the source of electrons due to neutral beam injection. In METIS, the parameter ${{f}_{n0}}$ is prescribed, one value for the limiter configuration and one for the divertor configuration; it can be read from external modelling made with a 2D edge code. It is worth to remark that the term $\frac{{\rm d}N}{{\rm d}t}$ accounts also for gas puffing and pumping as the density profile is controlled through the prescribed line-averaged density.

Finally, from the neutral source the ionisation heat sink is computed.

In the D–T mixture case, the Helium density is calculated as the sum of three terms, one representing the Helium resulting from the D–T fusion reactions $\left\langle {{n}_{{\rm He},{\rm ash}}} \right\rangle$, the second the fast alphas stored in the plasma $\left\langle {{n}_{{\rm He},{\rm fast}}} \right\rangle$ and the third representing the He coming from fuelling

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left\langle {{n}_{{\rm He},{\rm fuel}}} \right\rangle :\left\langle {{n}_{{\rm He}}} \right\rangle =\left\langle {{n}_{{\rm He},{\rm ash}}} \right\rangle +\left\langle {{n}_{{\rm He},{\rm fast}}} \right\rangle +\left\langle {{n}_{{\rm He},{\rm fuel}}} \right\rangle .\nonumber \end{align} \tag{ D.1 }$$

The first term ${{n}_{{\rm He},{\rm ash}}}$ is computed with a 0D equation for ashes accumulation:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{{\rm d}\left\langle {{n}_{{\rm He},{\rm ash}}} \right\rangle }{{\rm d}t}=~-\frac{\left\langle {{n}_{{\rm He},{\rm ash}}} \right\rangle }{{{\tau }_{{\rm He}}}}+~\frac{{{S}_{\alpha }}}{{{V}_{{\rm p}}}} \nonumber \end{align} \tag{ D.2 }$$

where ${{S}_{\alpha }}$ is the source of thermal alpha particles created by D–T reactions and ${{\tau }_{{\rm He}}}$ is the helium effective confinement time, provided either by a scaling law or by the following law adapted from [104]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\tau }_{{\rm He}}}=~{{f}_{{\rm He }~ }}{{\tau }_{E+~}}\frac{R}{1-R}{{\tau }_{{\rm ne}}}~\nonumber \end{align} \tag{ D.3 }$$

which allows to take into account the recycling of helium ashes.

The second term is computed from the fast alpha plasma energy content (that is an output of analytical Fokker–Planck solution used to compute the power deposition on electrons and ions):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left\langle {{n}_{{\rm He},{\rm fast}}} \right\rangle =~\frac{{{E}_{\sup ,\alpha }}}{2e3.56~{{10}^{6}}~{{V}_{{\rm p}}}~}.\nonumber \end{align} \tag{ D.4 }$$

The helium density profile is calculated assuming the same shape as the electron density profile, but with an edge value depending on ashes accumulation (more peaked profile and lower edge density due to a source located at the centre of the plasma):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{n}_{{\rm He}}}\left(t,x \right)=~{{C}_{{\rm He}}}\left(t \right)\left[ {{n}_{{\rm e}}}\left(t,x \right)-~{{n}_{{\rm e}}}\left(t,1 \right) \right]+~{{n}_{{\rm He},{\rm a}}}\left(t \right)\nonumber \end{align} \tag{ D.5 }$$

where ${{C}_{{\rm He}}}$ is set to have the desired $\left\langle {{n}_{{\rm He}}} \right\rangle$ value at each time slice and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{n}_{{\rm He},{\rm a}}}=~\left(\left\langle {{n}_{{\rm He},{\rm fuel}}} \right\rangle +~{{f}_{{\rm He},{\rm a}}}~\left\langle {{n}_{{\rm He},{\rm ash}}} \right\rangle \right)\frac{{{n}_{{\rm e}}}\left(1 \right)}{\left\langle {{n}_{{\rm e}}} \right\rangle } \nonumber \end{align} \tag{ D.6 }$$

where ${{f}_{{\rm He},{\rm a}}}\sim \frac{{{\tau }_{{\rm ne}}}}{{{\tau }_{{\rm He}}}}$

Remark: $\left\langle {{n}_{{\rm He},{\rm fast}}} \right\rangle$ have no contribution to edge density as fast alpha are confined to the centre of the plasma.

The second term $\left\langle {{n}_{{\rm He},{\rm fuel}}} \right\rangle$ is defined as a fraction of the average electron density. This fraction is prescribed by the user.

A specific treatment of the tungsten is included in METIS that users can switch on or off. This treatment is based on a model for the sputtering source at divertor target, a simple model for prompt-redeposition, a model for divertor screening and a simple model for tungsten accumulation in the plasma. Additionally, the source of electrons provided by local tungsten ionisation is taken into account and this effect is added to the electron density profile (that ensures to have the right electron density whatever is the tungsten accumulation). The temperature ${{T}_{{\rm e},{\rm t}}}$ and density ${{n}_{{\rm e},{\rm t}}}$ at the target and the temperature at LCFS must be provided in this case by the following model.

The profile of tungsten density in the plasma core is:

where ${{C}_{{\rm W}}}$, ${{C}_{{\rm W},{\rm offset}}}$, ${{\gamma }_{n~}}$, are user defined constants: ${{C}_{{\rm W}}}$ allows to take into account the strength of divertor source, ${{C}_{{\rm W},{\rm offset}}}$ allows to take into account the strength of other sources, ${{\gamma }_{n~}}$ the proportionality (1) or independence (0) of the tungsten density from electron density profile and is the strength of accumulation mechanism. This model for tungsten density assumes that each density at the LCFS is proportional to the flux [105].

The divertor response $G\left(t \right)$ is:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle G\left(t \right)={{Y}_{{\rm p}}}\left(t \right)~{{S}_{{\rm r}}}\left(t \right)~(1-~{{F}_{{\rm r}}}\left(t \right))\nonumber \end{align} \tag{ D.8 }$$

where ${{Y}_{{\rm p}}}\left(t \right)$ the sputtering factor, ${{S}_{{\rm r}}}\left(t \right)~$ the screening factor and ${{F}_{{\rm r}}}\left(t \right)$ the prompt-redeposition factor.

The sputtering factor reads:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{Y}_{{\rm p}}}\left(t \right)=~\underset{j\in \left\{{\rm species} \right\}}{\mathop \sum }\,\frac{{{n}_{j}}\left(t,1 \right)}{{{n}_{{\rm e},{\rm a}}}\left(t \right)}~{{y}_{j\to {\rm W}}}\left({{E}_{j}}\left(t \right) \right) \nonumber \end{align} \tag{ D.9 }$$

where the sputtering yield ${{y}_{j\to {\rm W}}}$ is either given by the model from [72] or by [106]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{E}_{j}}\left(t \right)=&\,{{T}_{{\rm e},{\rm t}}}\left(t \right)+\left(2+{{\Gamma }_{{\rm H}}} \right)*~{{T}_{{\rm i},{\rm t}}}\left(t \right)-~0.5~{{Z}_{j~}}{{T}_{{\rm e},{\rm t}}}\left(t \right)\nonumber \\ &\ln \left(\frac{2\pi ~{{m}_{{\rm e}}}}{{{m}_{j}}}~\left(1+~\frac{{{T}_{{\rm i},{\rm t}}}\left(t \right)}{{{T}_{{\rm e},{\rm t}}}\left(t \right)} \right) \right).\nonumber \end{align} \tag{ D.10 }$$

The two points model provides only $~{{T}_{{\rm e},{\rm t}}}$, so we fixed $\newcommand{\e}{{\rm e}} {{T}_{{\rm i},{\rm t}}}\left(t \right)={{\eta }_{{\rm t}}}{{T}_{{\rm e},{\rm t}}}\left(t \right)$. The factor ${{\Gamma }_{{\rm H}}}$ is 1 for the isothermal case and $\frac{5}{3}$ for the adiabatic case.

The screening model reads [105]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{S}_{{\rm r}}}\left(t \right)={{{\rm e}}^{-{{\left(\frac{{{T}_{{\rm leak}}}\left(t \right)}{{{T}_{{\rm e},{\rm t}}}\left(t \right)} \right)}^{2}}}}~\nonumber \end{align} \tag{ D.11 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{T}_{{\rm leak}}}\left(t \right)=~3.8~{{10}^{-9}}~\overline{{{Z}_{W}}}\left({{T}_{{\rm SOL},{\rm eff}}}\left(t \right) \right)~\sqrt{{{n}_{{\rm e},{\rm t}}}\left(t \right)\Delta s\left(t \right)}.\nonumber \end{align} \tag{ D.12 }$$

With $\overline{{{Z}_{W}}}$ is the averaged tungsten ion charge number [107] computed for temperature ${{T}_{{\rm SOL},{\rm eff}}}\left(t \right)={{f}_{{\rm SOL},{\rm ef\!f}~}}{{T}_{{\rm e},{\rm a}}}+$$\left(1-~{{f}_{{\rm SOL},{\rm eff}~}} \right)~{{T}_{{\rm e},{\rm t}}}$.

We assume that $\Delta s\left(t \right)$ is of the order of the neutral ionisation length near the divertor target:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Delta s\left(t \right)=~\frac{\sqrt{\frac{2e~{{T}_{{\rm e},{\rm t}}}\left(t \right)~}{{{m}_{{\rm eff}}}\left(t \right)}}}{{{n}_{{\rm e},{\rm t}}}\left(t \right)\left({{\left\langle \sigma v \right\rangle }_{{\rm i},{\rm e}}}\left(t \right)+{{\left\langle ~\sigma v \right\rangle }_{{\rm i},{\rm i}}}\left(t \right)+~{{\left\langle \sigma v \right\rangle }_{{\rm cx}}}\left(t \right) \right){{\mu }_{{\rm t}}}\left(t \right)}\nonumber \end{align} \tag{ D.13 }$$

where ${{\left\langle \sigma v \right\rangle }_{{\rm i},{\rm e}}}$, ${{\left\langle \sigma v \right\rangle }_{{\rm i},{\rm i}}}$, ${{\left\langle \sigma v \right\rangle }_{{\rm cx}}}$ are respectively electron ionisation, ion ionisation and charge exchange cross-sections for the main plasma species. ${{\mu }_{{\rm t}}}\left(t \right)$ is the angle of magnetic field line. The prompt-redeposition of tungsten ion can be switched on or off. If it is off then ${{F}_{{\rm r}}}\left(t \right)=~0$. Otherwise METIS uses a simplified model based on the ratio between tungsten ionised one time Larmor radius and ionization length [108] and a fit of redeposition computed by the ERO code results [109]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{F}_{{\rm r}}}~=~{{f}_{{\rm ERO}}}\left(\frac{{{l}_{W,{\rm ion}}}}{{{\rho }_{{{W}^{+}}}}} \right).\nonumber \end{align} \tag{ D.14 }$$

Here is the piecewise cubic Hermite Interpolated value of tabulated ERO results:

$x$	0	0.2	0.6	2.2	3.2	6	1000
${{f}_{{\rm ERO}\left(x \right)}}$	1	0.85	0.65	0.3	0.25	0.15	0.1

and assuming emergent tungsten ions have $~{{T}_{{\rm e},{\rm t}}}$ kinetic energy we have ${{l}_{{\rm W},{\rm ion}}}=~\frac{\sqrt{\frac{2e~{{T}_{{\rm e,t}}}}{{{m}_{{\rm W}}}}~}~}{{{\left\langle \sigma v \right\rangle }_{{\rm W,ionised}}}\left({{T}_{{\rm e,t}}} \right)~{{n}_{{\rm e,t}}}}$ and ${{\rho }_{{{{\rm W}}^{+}}}}$ is the tungsten Larmor ionised one time radius at divertor target.

The accumulation factor is chosen as the neoclassical expression in cylindrical geometry [40] with ad-hoc modification to take into account plasma rotation and heating decontamination effects:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{{\rm d}I}{{\rm d}x}={{C}_{\Omega ,{\rm W}}}~{{u}_{\Omega }}+~\frac{\frac{{\rm d}{{n}_{{\rm i}}}}{{\rm d}x}}{{{n}_{{\rm i }~ }}}+\left(1-\frac{1}{\overline{{{Z}_{{\rm W}}}}}-\frac{3}{2}-~{{C}_{{\rm H},{\rm W}}}~{{h}_{{\rm e,i}}} \right)\frac{\frac{{\rm d}{{T}_{{\rm i}}}}{{\rm d}x}}{{{T}_{{\rm i}}}} \nonumber \end{align} \tag{ D.15 }$$

where ${{C}_{\Omega ,{\rm W}}}$ and ${{C}_{{\rm H},{\rm W}}}$ are user-defined constants that drive the contamination and decontamination processes due to, respectively, toroidal plasma rotation and additional heating. The rotation term $~{{u}_{\Omega }}$ is [110, 111]:

$~{{u}_{\Omega }}=~-~\frac{\frac{{\rm d}\Omega }{{\rm d}x}}{{{R}_{{\rm axe}}}~{{v}_{{\rm th},{\rm W}}}}$ where ${{v}_{{\rm th},{\rm W}}}$ is the thermal velocity for tungsten.

This model is completed with a two points model (i.e. including neutral friction, radiative losses, momentum losses, kinetic correction and supersonic correction [105]) that allows to predict LCFS electron temperature and electron density and temperature at the outer divertor target knowing the heat flux crossing the LCFS and the electron density at LCFS. The scrape-off-layer (SoL) width is provided by a scaling law [96]. Eventually, this model takes into account counter-reaction due to tungsten source from sputtering on radiative fraction both in core plasma, in SoL and divertor region.