Testing validity of 1D models for impurity fraction scaling for divertor detachment with EDGE2D-EIRENE

Large tokamaks of near-reactor or reactor size, such as ITER [1] and DEMO [2], are expected to operate at high plasma densities with partially detached divertors in order to reduce power loads on divertor plates. The detachment is assisted by injecting (also 'puffing' or 'seeding') impurities (typically, light impurities such as N, Ne or Ar) radiating a large part of the power flowing to divertors. At the same time, penetration of these impurities into the main, or 'core' plasma (inside the separatrix), may present a problem for sustaining a high fusion reaction rate due to the dilution of ion species responsible for the fusion reaction (D⁺ and T⁺), as well as due to cooling of the core plasma, leading to a reduction of the fusion power. Optimal impurity concentrations required to exhaust power in the scrape-off layer (SOL) and divertor, while at the same time maintaining acceptable fusion power levels in the core, are therefore important parameters for a fusion reactor.

A number of analytical studies were recently carried out to develop scaling for impurity fractions (defined as the ratio of the total impurity density, including all ionization stages, over the electron density). The model of Kallenbach et al [3] allows the prediction of the seed impurity fraction required to achieve partial detachment for a given power flux into the SOL, P_SOL , in devices with a closed divertor similar to that of ASDEX Upgrade. The other three models use a simplified 1D geometry covering the SOL and divertor regions, assuming constant impurity fractions along magnetic field lines. They are partly [4] or fully [5, 6] based on dimensionless parameters characterizing an acceptable P_SOL (high enough for the total power input into the plasma, P_in , to exceed the minimum power level required for the plasma to stay in the high-energy confinement mode, or 'H-mode'), sufficiently high plasma density necessary for high fusion yield, and the 'radial' decay length of the parallel power flux in the SOL in the direction perpendicular to the magnetic surfaces. In particular, the Reinke model [5] solves a complete set of equations in order to eliminate an explicit dependence on P_SOL by replacing it with the dependence on the f_LH coefficient, which is the factor of excess of P_in over Martin's experimentally established H-mode power threshold [7]. This is justified by the link between c_Z and the impurity radiation power P_rad,imp (or P_rad,N in this paper due to the injection of nitrogen), which in detached conditions approaches the power flux into the SOL. The explicit dependence on the separatrix electron density n_e,sep is replaced by the dependence on the Greenwald density limit [8] fraction f_GW (as in [4]), here defined as n_e,sep /n_e,GW . The radial decay length of the parallel power flux in the SOL ${\lambda _q}$ is related to the experimentally established Eich scaling [9]. In a simple 1D SOL (and divertor) model, ${\lambda _q}$ represents the SOL width.

The most recent of the models listed above is the Huber–Chankin (HC) model [6]. Among other models, the Reinke model [5] is ideologically the closest to the HC model. Unlike [5], the [6] model contains a hydrogen isotope mass scaling (covering H, D and T isotopes) of the main working gas used in the discharge. Correspondingly, Martin's scaling for the H-mode power threshold P_LH is corrected by the inclusion of the dependence on the isotope mass (inversely proportional to it) first established experimentally in [10]. Also, motivated by tracing the impurity fraction dependence on the isotope mass, the HC model uses Goldston's heuristic drift-based (HD) model for ${\lambda _q}$ [11] instead of Eich's model, which does not have isotope mass dependence. The HD model for ${\lambda _q}$ is in a reasonable quantitative agreement both in absolute magnitude and in scaling with the experimental data [11]. Finally, the f_GW dependence in [5] is replaced by the dependence on the ratio of n_e,sep over the critical separatrix electron density n_e,crit following from Goldston's model for the density limit caused by the MHD instability in the SOL close to the separatrix [12]. Its use eliminates the need to relate the separatrix density, one of the key parameters in this study, to the Greenwald density limit, which depends on the line-averaged density. It has to be pointed out that in models [4–6] the only source of volumetric power loss is impurity radiation.

The present study is dedicated to testing the impurity fraction scaling following from the HC model using the EDGE2D-EIRENE code package combining the coupled plasma 2D edge fluid code EDGE2D with the neutral Monte-Carlo solver EIRENE [13–15]. The code package covers the outer core, SOL and divertor regions. Of particular interest is the scaling for the impurity density c_Z including all of its charge states, n_Z , divided by the electron density n_e and averaged over all cells of the first poloidal ring (of the EDGE2D grid) outside the separatrix, from X-point to X-point. In this averaging the sum of f_Z cell values multiplied by the poloidal cell sizes was divided by the sum of the poloidal sizes. As pointed out in the paper on the physics basis for the first ITER tungsten divertor, containing a detailed description of simulations with the edge 2D fluid code SOLPS coupled with EIRENE [16], this is the most convenient quantity for coupling SOL to a core model, and has been the key point of the SOLPS-4.3 model for ITER since 2003. This quantity is typically quoted in percentage terms.

As will be discussed later in the paper, assumptions used in the HC and similar models based on a simple 1D geometry may come into conflict with 2D distributions of plasma parameters in the SOL and divertor. Since no experimental data are used in this paper, the present study represents a model (1D) to code (2D) comparison. Regardless of the exact assumptions of a 1D model, certain features, such as the neglect of neutral ionization processes, the very different plasma parameters in the SOL and divertor, and the absence of friction and thermo-forces exerted on impurities by the main ions potentially leading to significant differences between impurity fractions in the SOL and divertor (while a spatially constant impurity fraction is assumed in a model), can render the model results unreliable. In this respect, the present study is, more generally, a test of the ability of simple 1D models to predict impurity fractions with a certain degree of reliability when confronted with the 2D effects of a real magnetic geometry, which are omitted in 1D models but included in much more detailed 2D models.

In this paper, section 2 describes the HC model and the strategy of its verification by the EDGE2D-EIRENE code. The setup of the EDGE2D-EIRENE ('code') cases is described in section 3. The predictions of the HC model are compared with the code results in section 4. The power balance in the code cases is analyzed in section 5. Consistency between the full HC model scaling for c_Z , including some of its parameters that cannot be predicted, and the code output is checked in section 6. The causes of the discrepancies between the model predictions and code results are analyzed in section 7. The results of this work are summarized in section 8.

The final scaling for c_Z derived in [6] (equation (37)) is replicated here without a numerical coefficient and with some re-arrangement of terms:

$$\begin{equation}{c_Z} \propto {{{B^{1.47}}} \over {{L_{INT}}}}{{q_c^{0.39}} \over {{A^{1.38}}}}{{f_{LH}^{1.38}} \over {{f_{ne}}{f_{sep}}}}{1 \over {G_A^{1.19}G_Z^{0.38}}} \times {{{R^{1.59}}{{\left(1 + {k^2}\right)}^{0.14}}} \over {{\varepsilon ^{0.08}}}}.\end{equation} \tag{ 1 }$$

Here, ${c_Z} = {n_Z}/{n_e}$, and the ratio is assumed to be constant along the field lines in the model. B is the toroidal field, ${L_{INT}} = \int_{{{\text{T}}_{e{\text{,}}t}}}^{{T_{e{\text{,}}u}}} {{L_Z}} \left( {T_e^{\prime}} \right)\sqrt {T_e^{\prime}} {\text{d}}T_e^{\prime}$, where L_Z is the cooling rate coefficient due to impurities, T_e,t and T_e,u are the target and upstream electron temperatures, q_c is the cylindrical safety factor, and A is the hydrogen isotope mass. f_LH = P_in /P_LH is the ratio of the input power by the H-mode power threshold, f_ne = n_e,sep /n_e,crit , with the expression for n_e,crit given by equation (28) of [6]; the scaling for n_e,crit relevant to the present study is given by equation (A4) in appendix . f_sep is the ratio of n_e,sep by the electron line-averaged density ${\bar n_e}$ figuring in Martin's H-mode power scaling (where we add the dependence on the isotope mass A to this scaling, since it was used in the derivation of equation (1), and remove the numerical coefficient):

$$\begin{equation}{P_{LH}} \propto \frac{1}{A}\bar n_e^{{\text{0}}{\text{.72}}}{B^{{\text{0}}{\text{.8}}}}{S^{{\text{0}}{\text{.94}}}} \propto \frac{{n_{e{\text{,}}sep}^{{\text{0}}{\text{.72}}}{B^{{\text{0}}{\text{.8}}}}}}{A}.\end{equation} \tag{ 2 }$$

In the above scaling, in its second proportionality, we replace ${\bar n_e}$ with n_e.sep and drop the ${f_{sep}} = {n_{e,sep}}/{\bar n_e}$ coefficient, implicitly assuming proportionality between these two quantities. This is a reasonable assumption since the power threshold P_LH is known to be dependent on the edge rather than the line-averaged density (see e.g. [17, 18]) and because we will be using the EDGE2D-EIRENE code to verify the scaling equation (1), and this code covers only the edge plasma. Hence, retaining the f_sep coefficient cannot potentially bring any additional information. Also, in the second proportionality of equation (2) we drop the dependence on the area S of the flux surface containing the magnetic separatrix since the same EDGE2D (JET) grid will be used in all EDGE2D-EIRENE cases considered here.

The G_A and G_Z parameters in equation (1) are related to the exact composition of the ion component:

$$\begin{equation}{G_A} = \frac{{2\bar A}}{{1 + \bar Z}},{G_Z} = \frac{5}{{{Z_{eff}} + 4}}\end{equation} \tag{ 3 }$$

where $\bar A$ is the average ion mass, including all ionization stages of impurities, $\bar Z$ is the average ion charge, again including all ionization stages of impurities, and Z_eff has its usual meaning.

Finally, the dependencies on the major radius R, plasma elongation k and inverse aspect ratio $\varepsilon = a/R$ in equation (1) are grouped together in terms separated by the multiplication sign '×'. They will not be used in the present study, which is based on EDGE2D-EIRENE runs using only one computational grid, as was pointed out above. The dependence on these parameters will therefore be omitted in the rest of the paper.

In order to use the EDGE2D-EIRENE output results for comparison with the scaling equation (1), in this work another parameter was introduced into the set of equations, and the dependence of c_Z on it was traced in the derivations. This parameter is the ratio of the heat flux decay length (or the SOL 'power width' used in the 1D model) to that following from the HD model: ${f_{HD}} = {\lambda _q}/\lambda _q^{HD}$, where ${\lambda _q}$ is the power decay length equal to 2/7 ${\lambda _{Te}}$, and $\lambda _q^{HD}$ is given by equation (6) of [11]. Since it is difficult to match ${\lambda _q}$ in the code calculations to $\lambda _q^{HD}$, the equations in [6] were re-derived to track the dependence on this additional parameter, yielding the dependence $f_{HD}^{ - 1.43}$ that has to be added to the right-hand side of equation (1); see equation (6).

Special attention should be paid to the G_A dependence in equation (1). For very low impurity concentrations, $\bar A \approx A$ and $\bar Z \approx 1$ (for hydrogen isotopes as the main gas), and the dependence of ${c_Z}$ on ${G_A}$ scales as ${G_A} \propto A$ against A, resulting in the scaling ${c_Z} \propto {A^{ - 1.38}}$ becoming much stronger: ${c_Z} \propto {A^{ - 1.38}}G_A^{ - 1.19} \propto {A^{ - 1.38}}{A^{ - 1.19}} = {A^{ - 2.57}}$. For the impurity content relevant to this study, ${G_A}$ has a weaker dependence on A. For example, for nitrogen as the only impurity and ${c_Z} = 0.03$ (3%), ${G_A} \propto {A^{ - 0.88}}$ is found, and this dependence is further weakened with the increase in the nitrogen impurity density. Nevertheless, as will be demonstrated later in the paper, the dependence of G_A on A in this study for most of the EDGE2D-EIRENE code cases is strong enough for it to be included in the main part of the c_Z scaling. We therefore include it by introducing the additional parameter dependence ${A^{ - 1.19}}$ plus a redefined correction factor

$$\begin{equation}{G_{AA}} = \frac{{2\bar A}}{{A\left( {1 + \bar Z} \right)}}.\end{equation} \tag{ 4 }$$

We now present the final c_Z scaling to be used in this paper, including the dependence on f_HD relevant to the present study:

$$\begin{equation}{c_Z} \propto {{{B^{1.47}}q_c^{0.39}} \over {{A^{2.57}}}}{{f_{LH}^{1.38}} \over {{f_{ne}}}} \times {f_{corr}},\end{equation} \tag{ 5 }$$

where

$$\begin{equation}{f_{corr}} = {1 \over {{L_{INT}}}}{1 \over {G_{AA}^{1.19}G_Z^{0.38}}}{1 \over {f_{HD}^{1.43}}}.\end{equation} \tag{ 6 }$$

Note that the dependence on f_sep does not figure in equations (5) and (6) since any reference to the line-averaged density was removed in the present work, as discussed above.

Dependencies on A, B, q_c, f_LH and f_ne will be considered to be the main dependencies. They can be used for predictions of c_Z values. The first three parameters, A, B and q_c , will be varied in EDGE2D-EIRENE runs. Other dependencies should be considered as corrections; see equation (6). They cannot be used at the stage of setting input parameters for EDGE2D-EIRENE runs and can only be extracted from the code output. For example, despite the HC model using $\lambda _q^{HD}$ as the power decay length, the exact match ${\lambda _q} = \lambda _q^{HD}$ is difficult to achieve in the code modelling; hence, the correction coefficient f_HD extracted from the output of the code cases was introduced in the HC model, as explained above. This applies to all f-coefficients. As for f_LH , it cannot have the same value for all cases, because the H-mode threshold power depends on the density, which, in turn, can deviate from the aimed dependence equation (A2), hence the appearance of the case-to-case f_ne variation. It has to be stressed that the actual n_e,sep value achieved in each code case cannot, strictly speaking, be 'set'; it is rather the 'target' value. The difference between these two values is given by the f_ne correction factor.

The starting EDGE2D-EIRENE case (also referred to as the 'code case') in this study is a simulated radiative partially detached inter-Edge Localized Modes (ELM) H-mode JET plasma with an ITER-like wall in a vertical divertor configuration (with both strike points on vertical targets) based on one of the cases in [19]. The main working gas was deuterium, which was puffed from the inner divertor above the strike point. To cool the plasma, nitrogen (N) was evenly puffed as molecules (N₂) from the outer, low magnetic field side, between the entrance to the outer divertor and the outer midplane (OMP). Beryllium (Be) sputtered from the wall was also part of the ion mix in the simulation; however, it had very little effect on the plasma dominated by D (deuterium, puffed as D₂ molecules) and N puffs. The plasma density was controlled by the D puff and pumps positioned in the corners of the EDGE2D grid in the divertor, with the constant pump albedo. For neutral atoms and molecules, and impurity atoms, the Kotov-2008 EIRENE model [20] was used. Further details of the experimental data and simulations can be found in [19].

The EDGE2D grid used for the simulations can be found in [21], figures 2 and 3. The catalogued no-drift case selected from this study, ajarvin/edge2d/jet/85274/feb2116/seq#3, corresponded to a partially detached plasma with 8 MW of input power into the grid and 5 MW of nitrogen radiation. In the present study, Be was removed from the plasma ion mix and the case was re-run under the latest EDGE2D version. The case was then pushed towards the maximum nitrogen radiation P_rad,N = 6.03 MW by increasing the N puff. The catalogued case is alexc/edge2d/jet/85274/jul2223/seq#1. This case is labelled as case D (for deuterium) here and is also often referred to as the reference D case.

The OMP profiles of the electron and ion temperatures, T_e and T_i , and electron density n_e for this case are shown in figure 1.

Zoom In Zoom Out Reset image size

The target profiles of power density P_target, T_e and n_e for the reference D case are shown in figure 2. Very low P_target and T_e values around the strike points indicate that most of the input power in this case was radiated. Details of the power balance for this case will be presented in section 5.

Zoom In Zoom Out Reset image size

For comparison with the isotope scaling of the HC model, two more EDGE2D-EIRENE code cases were run, with deuterium replaced by hydrogen (case H) and tritium (case T). In addition, based on the D case, for comparison with B and q_c scalings of the HC model, four more cases were run; with higher and lower toroidal field B (cases HB and LB, respectively), and with higher and lower q_c than in the D case (cases Hq and Lq, respectively). The input parameters for these cases are consistent with dependences on fixed dimensionless coefficients f_LH and f_ne = n_e,sep /n_e,crit . In appendix the scaling for the input power, equation (A3), and separatrix density, equation (A4), for setting up additional code cases are derived.

Variation of the input power P_in has a profound effect on the edge plasma. It was varied by a factor of 1.3, with P_in in the HB case being larger and in the LB case lower than in the D case. (Keep in mind that the HB, LB, Hq and Lq cases all had deuterium as the main working gas). According to the P_LH scaling in equation (2), the factor 1.3 variation in P_in (higher and lower by this factor) requires factor 1.22 variation of B in the code cases. The cylindrical safety factor q_c was varied by a factor of 1.5, again, higher and lower by this factor, in the cases Hq and Lq, respectively.

Cases with higher and lower q_c used the same EDGE2D grid, but the code plasma transport equations were changed. In EDGE2D, the effect of the higher q_c by a factor of 1.5 is achieved by specifying the extra namelist parameter BTHFACT, which, when set to 1.5, multiplies the ratio of poloidal to total magnetic fields by this factor in transport equations. The lower q_c is achieved, correspondingly, by setting BTHFACT = 0.667.

It has to be pointed out that the value of B in all code cases by itself plays no role since, in the EDGE2D-EIRENE model described here, drifts were switched off. The influence of the B and q_c variation on the code case setups, as well as the influence of the isotope mass A variation comes via dimensionless parameters. For example, using H instead of D increases P_in by a factor of 2 due to the inverse isotope mass scaling of the H-mode power threshold given by equation (2). Since the isotope mass A also figures in the scaling for the critical separatrix electron density n_e,crit (see equations in [6], as well as equations (A2) and (A4), where it is denoted as n_e,sep ), n_e,sep must be changed by varying the hydrogen isotope gas puff (GP) level, and similarly the power decay length $\lambda _q^{HD}$, which also has dependence on A.

The decay length ${\lambda _q}$ cannot be easily controlled in the code run. In this study, its variation was achieved by varying anomalous perpendicular transport coefficients for electron and ion power fluxes. According to the models in [4–6], the decay length of the parallel power flux in the SOL is attributed to the decay length of the electron conductive power flux given by ${\lambda _q} = 2/7 \times {\lambda _{{T_e}}}$, where ${\lambda _{{T_e}}}$ is the radial decay length of the electron temperature. In order to be able to change ${\lambda _q}$ in the code, the electron perpendicular heat conduction ${\chi _e}$ has to be varied according to the scaling equation (A9) derived in appendix. The factor of the variation in ${\chi _e}$, following from this scaling, was applied to each radial position of the EDGE2D grid. (In EDGE2D input panels, transport coefficients are specified for a few radial positions at the OMP. Linear interpolations onto cell center positions between each two successive transport coefficients are done using the code. Therefore, variation of the entered transport coefficient values by a certain factor ensures that for all cells of the EDGE2D grid, variations in transport coefficients are calculated using the same factor.) The same factor was also applied to the ion heat conduction coefficient ${\chi _{i \bot }}$, but this procedure was not extended to the particle diffusion coefficient ${D_ \bot }$ since, according to models [4–6], the n_e decay length has no impact on c_Z scaling, and only ${\lambda _{{T_e}}}$ is taken into account when defining the SOL (power) width.

Generally, using variable ${\chi _{_ \bot }}$ coefficients is not absolutely necessary, since any deviation of ${\lambda _q}$ from $\lambda _q^{HD}$ extracted from the code output can be included in the correction factor ${f_{HD}}$. Still, varying ${\chi _{_ \bot }}$ makes the code case parameters and profiles closer to those of the used model, reducing the correction factor f_HD .

In this section we compare predictions of the c_Z scaling given by equation (5), but without the f_corr coefficient, with c_Z values extracted from the code runs for plasma parameters along the first poloidal ring outside of the separatrix (ring s01, using EDGE2D nomenclature), from the outer to inner target plates in the divertor, in line with cell counting in EDGE2D. The radial profile effects relating the code output to the simple 1D model scaling are accounted for by including the radial power decay length ${\lambda _q}$ at the OMP position, corresponding to the SOL width of the HC model scaling.

The separatrix density correction factor f_ne is also included in this predictive part of the scaling since the electron density at the OMP position is the input parameter for analytical scaling, while in the code results some deviations from the 'setup' n_e values, following from equation (A4), occurred.

The input power into the EDGE2D grid P_in in setting up the EDGE2D-EIRENE cases follows equation (A1), which, after elimination of the n_e,sep dependence, yields the scaling (A3) (which assumes that n_e,sep follows the scaling equation (A4)). Since the separatrix electron density in the code cases did not match this scaling, the correction coefficient f_ne in equation (5) should also impact the f_LH coefficient via the dependence of P_LH on the separatrix density according to equation (2): ${P_{LH}} \propto n_{e,sep}^{0.72}{B^{0.8}}{S^{0.94}}/A$ (see the discussion on the density power threshold dependence in section 2). Hence, the density correction factor f_ne should also impact f_LH according to ${f_{LH}} \propto 1/n_{e,sep}^{0.72}$, since e.g. for higher code separatrix electron density compared to that following from the scaling equation (A4), the input power becomes lower compared to that used in the scaling for P_SOL ; hence, ${f_{LH}} \propto 1/f_{ne}^{0.72}$, giving the following total dependence on this parameter (and remembering that c_Z scales with $f_{LH}^{1.38}$):

$$\begin{equation}{c_Z} \propto {{{B^{1.47}}q_c^{0.39}} \over {{A^{2.57}}}}{1 \over {f_{ne}^{1.99}}} \times {f_{corr}}.\end{equation} \tag{ 7 }$$

The above scaling without f_corr will be used in this section to calculate the HC model predictive c_Z values instead of equation (5).

Table 1 contains relevant information on the EDGE2D-EIRENE cases. The catalogued case names are given by their suffixes, which must have the same prefix 'alexc/edge2d/jet/85 748/jul'. P_in is the input powers into the computational grid. B(T) is a toroidal magnetic field on axis, and q₉₅ is the safety factor at 95% of magnetic flux used as a proxy for q_c. T_e,sep and n_e,sep are the electron temperature and density at the OMP position. f_ne is the correction factor calculated by dividing code n_e,sep by the 'target' electron density following from the density scaling in equation (A4). Since the D case is considered to be a reference case on which other case parameters are often normalized, f_ne for this case is set to be equal to 1. '$c_Z^ *$ HC' is the nitrogen impurity fraction calculated according to the HC model (see equation (7) without the f_corr factor) and '$c_Z^ *$ code' is the same impurity fraction extracted from the EDGE2D-EIRENE solutions. An asterisk in '$c_Z^ *$' indicates that these values were normalized by those in the reference D case. It has to be acknowledged that such a normalization reflects the authors' assessment of the code's poor ability to predict absolute c_Z values owing to the number of simplifications adopted in the model, mostly related to its simple 1D geometry.

Table 1. Main parameters of EDGE2D-EIRENE code cases and comparison between the HC model predictions and code results.

Case label	Case name (suffix)	Pin MW	B T	q95	Te,sep (eV)	ne,sep (109m−3)	fne	$c_Z^ *$HC	$c_Z^ *$code
D	2223/seq#1	8.0	2.67	3.20	89.9	89.9	1	1	1
H	1823/seq#1	15.1	2.67	3.20	140.7	140.7	1.56	2.45	2.73
T	1923/seq#4	5.52	2.67	3.20	70.6	70.6	0.94	0.40	0.24
HB	1923/seq#1	10.4	3.26	3.20	100.5	100.5	1.13	1.05	1.18
LB	1923/seq#2	6.16	2.19	3.20	64.6	64.6	0.98	0.78	1.36
Hq	2023/seq#1	5.90	2.67	4.80	64.2	64.2	0.94	1.32	1.34
Lq	1923/seq#3	10.8	2.67	2.13	91.5	91.5	1.01	0.84	1.79

The absolute values for (poloidally averaged along the separatrix, from X-point to X-point, along ring s01) code ${c_Z}$ values in % are (from D to Lq case): 4.47, 12.21, 1.06, 5.28, 6.08, 6.01, 8.00.

As can be seen from this table, apart from the high correlation between the model and code $c_Z^ *$ values for the D, H and T cases, there is no correlation (or rather, a slightly negative correlation of −0.34, to be precise) between these values for cases in deuterium plasmas (cases D, HB, LB, Hq and Lq). In particular, between the pairs of cases HB, LB and Hq, Lq, trends of the $c_Z^ *$ values are opposite. There is also a large difference between the $c_Z^ *$ values from the HC model in cases LB and Hq, while in the code the $c_Z^ *$ values are very close to each other.

The good correlation between the model and code $c_Z^ *$ values in the isotope scan (H, D and T) can be attributed to the very large difference in P_in in these cases caused by the use of the P_LH ∼ 1/m_i dependence. The P_in variation among the isotope cases is significantly larger than that among the five cases in deuterium. Since in the modelling all cases were pushed to the maximum radiated power, mostly on N, these cases had the largest difference in N radiation (see P_rad,N values in table 2 in the next section). The large P_rad,N difference between the code cases does not necessarily guarantee a large c_Z difference, but at least it is a strong argument to consider this to be a plausible explanation for the above-mentioned correlation.

Table 2. Parameters related to the power balance in the code cases.

Case label	PSOL MW	Prad,N (MW), Prad,N /PSOL	Pother MW	Ptar, Pwall MW	$\tilde f_{LH}^*$ (Prad,N )
D	7.84	5.42, 69%	1.13	0.76, 0.56	1
H	14.11	10.71, 76%	0.77	2.60, 0.14	1.21
T	5.46	2.87, 53%	1.13	0.59, 0.89	0.89
HB	8.74	5.16, 59%	2.19	1.12, 0.41	0.85
LB	5.39	4.44, 82%	0.52	1.10, 0.32	1.23
Hq	4.07	2.59, 64%	0.52	0.12, 0.85	0.73
Lq	8.6	6.28, 73%	1.73	0.56, 0.09	1.34

Two cases in table 1 stand out. The first is the reference D case. One would expect its code $c_Z^ *$ value to be between the HB and LB case values, as well as between the Hq and Lq case values, since the B and q_c values for the D case are between those used in the above-mentioned pairs of cases. At least, this is expected from the HC model. In reality, the code $c_Z^ *$ value in this case is lower than in all of the other above-mentioned deuterium cases. Significant efforts were applied to raise ${c_Z}$ in this case by running numerous cases with different scenarios regarding the application of D and N GPs. None of them was successful, and the code ${c_Z}$ value remained below those in all four other deuterium cases. This may be related to effects not included in the simple 1D HC model, as will be further discussed in section 7.

The second unusual case is the H case, as its density correction factor f_ne (=1.56) strongly exceeds unity. Despite efforts made to reduce n_e,sep by gradual case-by-case increases of the input power, D and N puff levels, this case could not be driven to the domain of plasma parameters that would satisfy the conditions of f_ne values ∼1 together with the large radiation fraction in the power balance. Without allowing a significant increase in n_e,sep it was not possible to avoid a 'burnthrough' of high-temperature plasma from upstream to the outer target, resulting in a significant fraction of the target power deposition in the power balance. This is obviously strongly at odds with the assumptions of the HC model, where the power flux through the separatrix is assumed to be radiated on nitrogen. Apparently, the input power of 15.1 MW was too high for the H case to be mostly radiated on N without allowing a significant increase in n_e,sep . As will become clear below in section 6 (see table 3), the H case plasma had a very high content of N ions of various ionization stages. The ionization of nitrogen neutrals and ions creates too many electrons, leading to the uncontrolled n_e,sep rise.

Table 3. Output code parameters figuring in the full HC model scaling for detailed comparison between the model and the code.

Case label	${\lambda _{{T_e}}}$mm	$\lambda _q^ *$	${\bar A}$	${\bar Z}$	GAA	${Z_{eff}}$	${G_Z}$	$\lambda _ {q,HD} ^ *$	$f_ {HD}^{\,*}$
D	8.2	1	2.58	1.22	1.16	2.01	0.83	1	1
H	6.6	0.81	4.52	2.25	2.78	4.22	0.61	1.18	0.69
T	10.7	1.31	3.08	1.03	1.01	1.17	0.97	1.09	1.20
HB	7.7	0.94	2.74	1.28	1.20	2.23	0.77	0.88	1.07
LB	8.6	1.05	2.99	1.33	1.28	2.31	0.79	1.19	0.88
Hq	13.2	1.61	2.87	1.30	1.25	2.26	0.80	1.57	1.03
Lq	5.7	0.69	3.86	1.79	1.38	3.77	0.64	0.68	1.01

Failure to achieve the expected positive correlation of c_Z values between the model predictions and EDGE2D results may indicate that contributions from the quantities contained in the ${f_{corr}}$ term of equation (6), initially considered to be of secondary importance, may have a significant impact on the code predicted c_Z values. At the same time, in the code there exist several power loss channels, and the radiated power on N impurities is only one of them. In the next two sections, an analysis of power loss channels in the EDGE2D cases is performed, and the relative contribution of the impurity radiation is assessed, singling out this particular power loss and testing whether the model can reasonably describe c_Z code values when in the analysis the total power through the separatrix is replaced with the radiated power on N impurities (section 5). This analysis, together with the inclusion of the ${f_{corr}}$ multiplier in equation (6), is carried out in the following section 6.

Before moving to extraction from the code cases and/or calculations based on them of a number of parameters figuring in the f_corr factor equation (6), with the aim of achieving a better match between ${c_Z}$ values in the code and in the HC model, we would like to present some information regarding the power balance in the code cases in table 2. This is important since the code cases have a number of power loss channels that are not included in the HC model; for example, power losses on D (referred to as 'hydrogen' in EDGE2D) radiation or losses to the targets and the wall. If these power losses were to dominate the power balance, any comparison between the code and the (simple 1D) model would not make sense.

P_SOL in table 2 is the power flux actually delivered to the SOL and divertor regions from the core through the separatrix. P_rad,N is the nitrogen radiated power in the SOL and divertor. It, together with other power channels, can be extracted from the 'print' file created after the completion of each EDGE2D case, for each 'macro-zone' (SOL, inner and outer divertors connected to SOL along field lines, etc). It is this power that, strictly speaking, should be considered as P_SOL in the framework of the HC model's equations, since all power losses in this model are attributed to the radiation of impurities. As one can see from the third column of the table, P_rad,N accounts for more than 50% of the total power input through the separatrix; hence, the HC model, although not being detailed enough to cover all power losses, should at least be approximately applicable to the code case results.

P_other in the print files contains other volumetric power losses: atomic processes (neutral ionization and energy carried by neutrals), excited molecular states and dissociation energy, charge exchange energy, and 'hydrogen' (in our case deuterium) radiation. The remaining power losses are power fluxes to the target plates and the wall. None of the above power losses except for P_rad,N are included in the HC model.

$\tilde f_{LH}^*$ is the normalized (on its value in the D case) coefficient ${\tilde f_{LH}}$, defined as the ratio of P_rad,N inside the SOL and divertor over P_in in table 1. The P_in values in all cases ensure the same P_LH when taking into account only A, B and q_c variations, but ignoring the actual n_e,sep in code cases that may differ from scaling equation (A4) by the correction factor f_ne . Since in this and the next sections we use P_rad,N inside the SOL and divertor instead of P_in as an input power, ${\tilde f_{LH}}$ is proportional to the ratio of P_rad,N over P_in :

$$\begin{equation}{\tilde f_{LH}} \propto \frac{{{P_{rad,N}}}}{{{P_{in}}}}.\end{equation} \tag{ 8 }$$

The dependence of c_Z on ${\tilde f_{LH}}$ (or $\tilde f_{LH}^*$) must be considered as part of the main dependence in the c_Z scaling, since, in the framework of the HC model, here we assume P_rad,N = P_SOL , where P_rad,N variation over cases can strongly deviate from that of P_in . The dependence on the correction factor f_ne was already included in the main part of the c_Z scaling; see equation (7). The c_Z scaling one should use when assuming P_rad,N = P_SOL can therefore be expressed as

$$\begin{equation}{c_Z} \propto \frac{{{B^{{\text{1}}{\text{.47}}}}q_c^{{\text{0}}{\text{.39}}}}}{{{A^{{\text{2}}{\text{.57}}}}}}\frac{{\tilde f_{LH}^{1.{\text{38}}}}}{{f_{ne}^{{\text{1}}{\text{.99}}}}} \times {f_{corr}},\end{equation} \tag{ 9 }$$

where f_corr is given by equation (6) as before.

We now proceed with the calculation of parameters figuring in the expression for f_corr in equation (6), which are entered in table 3. Parameters figuring in f_corr , as was pointed out earlier, cannot be used for the prediction of c_Z values following from the HC model since their calculation requires more output from EDGE2D-EIRENE runs. It is of interest, however, to check whether their inclusion in the c_Z scaling may improve its match with the code output. As with all other earlier calculated parameters, those figuring in table 3 are calculated at the OMP position.

${\lambda _{{T_e}}}$ is calculated assuming exponential T_e decay at the OMP position between poloidal rings s01 and s05 separated by 4.4 mm. ${\lambda _{q,HD}}$ is given by equation (A7), G_AA by equation (4) and G_Z by equation (3). G_AA and G_Z require calculation of $\bar A$, $\bar Z$ and ${Z_{eff}}$ dependent on the composition of the ion component, and both the hydrogen isotope and N ions, including their charges and masses.

The average ion mass $\bar A$ in table 3 is highest in the H case, where one would normally expect the lowest value due to the low mass of hydrogen. However, the very high impurity (N) content in this case, as was pointed out in section 4, makes the $\bar A$ value highest in this case. Also, the average ion charge ${\bar Z}$ and ${Z_{eff}}$ are highest in the H case. Using the values of these parameters, the calculation of G_AA and ${G_Z}$ is straightforward. Together with A, B, q_c and n_e,sep from table 1, they allow one to calculate ${\lambda _{q,HD}}$ and then ${f_{HD}}$, which is entered in table 3 in its normalized form as $f_{HD}^ *$.

The last parameter required for the calculation of f_corr is ${L_{INT}}$, which will be entered in table 4 in its normalized form as $L_{INT}^ *$. The integral for calculating L_INT was presented in section 2. Here, this integral was taken from zero up to the highest T_e along the poloidal ring s01, which is quite close to T_e at the OMP positions in table 1 (T_e,sep ). The radiative cooling rate coefficient L_Z figuring in the integral as a function of T_e was taken from [22] for the case of coronal equilibrium. This publication contains L_Z values for a number of impurities, including nitrogen, taken from the Atomic Data and Analysis Structure database. Figure 5 of [22] shows the L_z (T_e ) curve for nitrogen, with table 9 of this reference providing the polynomial fit coefficients. As can be seen in table 4, $L_{INT}^ *$ shows quite moderate variation across code cases, with the maximum deviation from the D case being higher by a factor of 1.22 for the H case. For deuterium cases HB, LB, Hq and Lq, the largest difference from the D case is for the LB and Hq cases; for them, $L_{INT}^ *$ is lower than in the D case by a factor of 1.15.

Table 4. The impact of the radiative cooling parameter L_INT and other parameters, considered to be of secondary importance and collected in the correction coefficient f_corr , on the comparison between c_Z values in the HC model and the EDGE2D-EIRENE code.

Case label	$L_{INT}^ *$	$c_Z^ *$ HC without fcorr	$f_{corr}^ *$	$c_Z^ *$ HC with fcorr	$c_Z^ *$ code
D	1	1	1	1	1
H	1.22	2.45	0.56	1.37	2.73
T	0.90	0.40	0.95	0.38	0.24
HB	1.06	1.05	0.87	0.97	1.18
LB	0.87	0.78	1.24	0.69	1.36
Hq	0.87	1.32	1.04	1.37	1.34
Lq	1.01	0.84	0.8	0.74	1.79

Table 4 also presents the model $c_Z^ *$ values copied from table 1 and labelled as '$c_Z^ *$ HC without f_corr ', for easy comparison with other similar quantities, and the f_corr coefficient following from the dependence in equation (6). It is clear that, apart from the non-deuterium cases H and T, f_corr , similarly to L_INT , does not exhibit significant variations compared to the D case, varying by a maximum factor of 1.14, being lower in the LB case.

Multiplication of '$c_Z^ *$ HC without f_corr ' values by f_corr gives the full $c_Z^ *$ values, calculated according to equation (9) with f_corr calculated according to equation (5), and labelled '$c_Z^ *$ HC with f_corr ' in the table. For easy comparison with the code, the values of $c_Z^ *$, labelled '$c_Z^ *$ code', were copied into table 4 from table 1.

Summarizing the results of this section, it is clear that using even the full c_Z scaling derived in model [6] does not provide a positive correlation (in fact, a small negative correlation of −0.38, similar to that when using model c_Z values without f_corr ) between the model and code c_Z for cases with deuterium as the main working gas. As for the subset of cases with different hydrogen isotopes (D, H and T cases), they exhibit at least a qualitative correlation between the model and code c_Z values.

In all code cases, a large difference in the local (inside each cell) ion impurity density fractions f_Z = n_N /n_e between the main SOL and divertor regions can be seen (in this paper the notation c_Z is reserved for the poloidal average of f_Z in the main SOL region). The poloidal profiles of the nitrogen ion density fraction f_Z in cases D, HB and LB are presented in figure 3. Despite having a small difference in $c_Z^ *$ of just a factor of 1.15 between the HB and LB cases, they have the largest difference in the average impurity fraction between the divertor and main SOL plasmas (see below) among deuterium cases. The impurity fractions are plotted against cell numbers from the outer to inner target instead of against the poloidal distance for the sake of better spatial resolution. The profiles in figure 3 show large drops in the f_Z values near the entrances to the inner and outer divertors. Such a strong difference between the f_Z values in the main SOL and in the divertor could be caused by the very different composition of impurity ions by charge states, with lower ionization states N^+Z dominating in the divertor. According to the EDGE2D output, the dominant contribution to f_Z in the divertor comes from singly ionized N ions, while in the main SOL N⁺ charge states 5–7 dominate.

Zoom In Zoom Out Reset image size

The strong difference between plasma parameters in the main SOL and divertor can be seen in the parallel (along B) profiles of T_e and n_e in figure 4 for the same cases as presented in figure 3. In these cases, which were aimed at achieving the maximum impurity radiation power, steep jumps in the T_e and n_e values across the boundary between the main SOL and divertors can be seen. The divertor plasma is, as usual, much 'colder' and 'denser' owing to the high neutral recycling at the target plates. The steepness of the transition is exacerbated by large cell volumes on both sides of the boundary, having a common node at the X-point position. For the LB case, which has one of the (two) lowest input powers among the deuterium cases, the transition is shifted towards the main SOL at the inner (high B) side by one cell.

Zoom In Zoom Out Reset image size

The total ion (nitrogen) impurity content along ring s01 in the divertor is calculated by multiplying the nitrogen impurity density n_N by cell volumes and summing up all cell contributions inside the inner and outer divertor regions to give N_N,div . The same is done for the electron density, giving the total number of electrons in the divertor along ring s01, N_e,div . The ratio f_Z,div = N_N,div /N_e,div gives the average ion impurity fraction in the divertor. The same is done for the main SOL region, giving f_Z,SOL = N_N,SOL /N_e,SOL . The ratio R_div/SOL

$$\begin{align}{R_{div/SOL}} = {{{f_{Z,div}}} \over {f{}_{Z,SOL}^{}}}\end{align} \tag{ 10 }$$

will be used to characterize the difference between the average impurity fractions in the divertor and main SOL. Its values are given in table 5 (left column, for ring s01).

Table 5.  c_Z values and divertor to main SOL ratios of impurity fractions.

Case label	${R_{div/SOL}}$ ring s01	${R_{div/SOL}}$ ring s02	${R_{div/SOL,grid}}$	${c_{Z}},\%$ ring s01	${c_Z},\%$ ring s02
D	7.91	6.95	4.29	4.47	4.95
H	2.53	2.32	1.44	12.21	11.9
T	18.2	15.1	6.50	1.06	1.32
HB	1.93	2.99	1.92	5.28	4.45
LB	4.15	5.15	2.78	6.08	7.51
Hq	2.38	5.02	4.38	6.01	5.80
Lq	2.10	2.66	3.19	8.00	7.65

Apart from the ${R_{div/SOL}}$ values being significantly above unity (with the average value for ring s01 among deuterium cases being 3.7), they also exhibit large case-to-case variations, ranging from 1.93 in the HB case to 7.91 in the D case, among deuterium cases, giving a factor of up to 4.1 variation. This is a significantly larger variation than that in the code c_Z (or $c_Z^ *$) values, with a maximum difference between them, when comparing D and Lq cases, of factor 1.79 (see table 4, right column). The higher case-to-case variation of impurity fraction ratios ${R_{div/SOL}}$ among deuterium cases, compared to the variation of impurity fractions c_Z , further indicates a high degree of decoupling between the main SOL and divertor regions regarding impurity content.

Ratios ${R_{div/SOL}}$ > 1 should contribute to the shift of the radiation power towards the divertor, in addition to the plasma parameters being more favourable for radiating low Z impurities there (higher density, lower temperature) compared to those in the main SOL. In the EDGE2D cases, the calculated nitrogen radiation is found to come predominantly from the divertor (summing up contributions from both the inner and outer divertors) plus two large volume cells in the main SOL with one of their nodes coinciding with the X-point. For cases D to Lq, as listed in the left column of table 5 for example, the fraction of the N radiation coming from this region, when divided into the total N radiated power from the whole s01 ring, is (in %): 98.6, 97.9, 99.0, 98.9, 87.3, 53.4, 98.5. Note the relatively low fraction of 53.4% in the Hq case compared to other cases, which must be related to the expansion of the high-density, low-temperature plasma from the divertor onto the main SOL cell adjacent to the X-point.

Despite the use of only one ring s01 parameters for comparison with the HC model being in agreement with the logic of the 1D HC model, as it yields the important (for the core plasma) c_Z quantity at the separatrix in the main SOL, it is of interest to check whether results based on this ring can be applied, at least qualitatively, to a wider region outside the separatrix. One might for example argue that, regarding the differences in impurity fractions between the main SOL and divertor, it is more relevant to use some more outward ring with cells positioned at a distance of ∼${\lambda _q}$ from the separatrix, better representing the middle of the 1D 'grid', since its width is determined by ${\lambda _q}$ in the HC model. The average ${\lambda _q}$ for all cases is 2.5 mm. The position of the center of the OMP cell belonging to ring s02 is 2.86 mm outside of the separatrix. This ring may therefore be regarded as an alternative to ring s01 in representing the impurity density distribution between the main SOL and divertor. The ${R_{div/SOL}}$ values for ring s02 are shown in table 5. They exhibit a very similar case-to-case variation as for ring s01. The correlation of ${R_{div/SOL}}$ ratios between these two rings is 98.0%.

Calculation of the ${R_{div/SOL}}$ ratio covering the whole grid (SOL width at the OMP position is 23.2 mm) requires multiplication of the total nitrogen ion density, as well as the electron density, by cell volumes inside of the divertor region, with the consequent division of the sum of the products over these cells by the same quantity calculated in the main SOL, yielding the ratio ${R_{div/SOL,grid}}$ (for the whole EDGE2D grid). These ratios are also presented in table 5. The correlation between these values and those for ring s01 is 83.3%, with the average value being 3.50 for the whole grid, significantly below the 5.70 value for ring s01. Still, this indicates that the impurity ion distribution between the divertor and main SOL can to some extent be applied to the whole divertor and main SOL regions.

Finally, table 5 also presents c_Z values calculated for ring s02, together with those for s01. The correlation between the two sets of data, covering all cases, is 97.6%.

A reasonably good match between the predictions of the HC model and code (EDGE2D-EIRENE) results on the poloidally averaged impurity fraction c_Z along the separatrix in the main SOL for D ('reference D case'), H and T cases has been obtained. The match can be explained by the inverse dependence of the H-mode power threshold on the hydrogen isotope mass of the main working gas, which in the code cases was set to be proportional to the input power into the EDGE2D computational grid, P_in . The difference in P_in among isotope (H, D, T) cases is significantly larger than among the five deuterium cases. For a higher P_in , a higher impurity density, achieved by higher nitrogen puffing, is required to radiate the input power, which is a precondition for detachment. A higher impurity content is associated with a larger impurity fraction ${c_Z}$. This may explain the positive correlation between ${c_Z}$ values in the D, H and T cases and the predictions of the HC model.

At the same time, no correlation was found between the model and code ${c_Z}$ values among the deuterium cases: the (reference) D case, cases with higher and lower toroidal magnetic field (HB and LB cases, respectively), and cases with higher and lower cylindrical q (Hq and Lq cases, respectively). The code cases exhibited larger, by a factor of 3.7 on average (using ${c_Z}$ values in deuterium cases), impurity fractions in the divertor than in the main SOL along poloidal ring s01. With the dominant part of the nitrogen impurity radiation coming from the divertor and two cells in the main SOL adjacent to the X-point, the model, which uses the assumption of a spatially (along field lines) constant impurity fraction, defined as the ratio of the sum of impurity species' densities over the electron density, and aiming to predict c_Z in the main SOL, becomes inadequate to describe the plasma parameters in the SOL and divertor. In addition to the large average ratio of the divertor to the main SOL impurity fractions, the code also gives a large (factor 4.1) case-by-case variation of this ratio (taking deuterium cases with the lowest and highest ratios).

The ultimate cause of the discrepancies between the model and code c_Z values must be the effects of the real magnetic equilibrium consisting of the main SOL and divertor regions, which have quite different plasma parameters created mostly by high neutral recycling in the divertor, plus the totally unjustified assumption that f_Z = const along field lines. The HC model and other similar models assume the simplest possible geometry and the absence of neutrals. Neutrals, apart from cooling the divertor plasma and raising the electron density by their ionization, create parallel plasma flows in both the main SOL and divertor, which is a 2D effect (see e.g. [23] chapter 15: 'Ionization-driven flow reversal', p 471). The friction force exerted on impurities by the main ions, together with the thermo-force caused by parallel ion temperature gradients, are known to strongly influence impurity density distributions along field lines.

Going back to the comparison of D, H and T cases, it is worth pointing out that the positive correlation between c_Z fractions in the model and the code could be achieved owing to exceptionally large differences in ${c_Z}$ values following from the code, ranging by a factor of 11.5 when comparing the H and T cases. This variation is even larger than the factor 7.2 variation in the ${R_{div/SOL}}$ ratio (again comparing the H and T cases), characterizing the difference between impurity fractions in the main SOL and divertor, which could be considered as 'scatter' influencing the c_Z values. Hence, the ${R_{div/SOL}}$ variation may not have dramatically affected the positive correlation of c_Z values mentioned above. Still, in terms of the predicted power of the HC model's scaling for the H, D and T plasmas (and leaving aside the usefulness of such a prediction, since tokamak reactors considered presently are supposed to operate with only one, approximately 50/50% D–T mixture), it is still unacceptable because of the strong difference in impurity fractions between the highly radiating divertor plasma, to which the HC and other 1D models largely apply, and the main SOL plasma at the separatrix position, for which 1D models aim to predict ${c_Z}$.

Retrospectively, conclusions on the origin of the large mismatches between the predictions of the 1D model and the results of the 2D code EGDE2D-EIRENE seem obvious, and may not have warranted such a detailed analysis as was performed in this work. However, the situation could have been different if, suppose, EGDE2D-EIRENE had incidentally produced an acceptable correlation between the model and code ${c_Z}$ values, since the amount of data for the statistical analysis is very limited. The detailed analysis performed here therefore seems justified since it focusses attention on the fact that 1D models to a large extent describe high-radiation plasmas in the divertor and close to the X-point in the main SOL, while what is required are predictions for ${c_Z}$ at the separatrix in the main SOL far away from the X-point, and there is a large degree of decoupling between impurity fractions in the two regions, judging by the large variation of ${R_{div/SOL}}$ ratios obtained in the code despite them being connected by power conduction from the main SOL to the divertor. In addition, important physics related to the poloidal distributions of impurity ions influenced by the balance between the friction and thermal forces between the main ions and impurities are absent in 1D models, and are replaced by a simplistic assumption of constant local impurity fraction f_Z along field lines. The 1D models could possibly be improved by accounting for the variation of local impurity fractions f_Z due to electron temperature variation along field lines via the different composition of impurity ions by their ionization states. This, however, would require also introducing the density dependence of the impurity composition (at present, the electron density does not explicitly figure in the 1D models). However, to account for parallel plasma flows and the friction force on impurities would require the introduction of neutrals in the models, ultimately pushing them towards becoming 2D models.

This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No. 101052200—EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.

All data that support the findings of this study are included within the article (and any supplementary files).

In this section we derive scalings for setting up the input parameters for the EDGE2D-EIRENE cases. The HC model's equations are derived using dimensionless parameters, such as the ratio of the input power by the H-mode power threshold and the ratio of the separatrix electron density by the critical electron density n_e,crit following from Goldston's model for the density limit caused by the MHD instability in the SOL [12]. The code cases' input parameters should, ideally, keep the dimensionless parameters unchanged, so that correction factors 'f' would not be too different from unities.

We start with the scaling for Martin's H-mode power threshold corrected for the isotope dependence, equation (2), leaving only parameters relevant for the case setups and ignoring the difference between the line-averaged and separatrix electron densities, as discussed earlier:

$$\begin{align}{P_{SOL}} \propto {{{B^{0.8}}} \over A}n_{e,sep}^{0.72}.\end{align} \tag{ A1 }$$

Equation (28) of [6] provides the scaling for the separatrix electron density:

$$\begin{align}{n_{e,sep}} \propto {{{B^{7/8}}P_{SOL}^{ - 1/8}} \over {q_c^{9/8}}}.\end{align} \tag{ A2 }$$

Substituting the n_e,sep scaling into equation (A1), one obtains

$$\begin{align}{P_{SOL}} \propto {1 \over {{A^{0.92}}}}{{{B^{1.32}}} \over {q_c^{0.75}}}.\end{align} \tag{ A3 }$$

Substituting the P_SOL scaling equation (A3) back into equation (A2) yields

$$\begin{equation}{n_{e,sep}} \propto {A^{0.12}}{{{B^{0.71}}} \over {q_c^{1.03}}}.\end{equation} \tag{ A4 }$$

As pointed out above, this separatrix density can be reached by varying the D GP. At the same time, since the ultimate goal is to reach the maximum radiated power in the code cases, the exact match of n_e,sep by varying both deuterium and nitrogen puffs may be difficult when problems such as the failure of code cases or the appearance of oscillations of plasma parameters are encountered. In such cases correction factors ${f_{ne}}$ normalized to the reference case D, are used. They account for the deviation of code separatrix electron density from the scaling (A4).

In the derivations of [6], the power decay length $\lambda _q^{HD}$, following from the HD model, was used. This decay length is given by equation (35) of [6]:

$$\begin{align}\lambda _q^{HD} \propto {\left( {{P_{SOL}}{\text{l}}_{{\text{||}}}^{\text{*}}} \right)^{{\text{1/8}}}}\frac{{q_c^{{\text{9/8}}}}}{{{B^{{\text{7/8}}}}}}{\left( {\frac{{2\bar A}}{{1 + \bar Z}}} \right)^{{\text{7/16}}}}{\left( {\frac{{{Z_{eff}} + {\text{4}}}}{{\text{5}}}} \right)^{{\text{1/8}}}}.\end{align} \tag{ A5 }$$

The dimensionless variable $l_{||}^*$ above is defined as $l_{||}^* = {L_{||}}/\pi {q_c}R$, while the parallel connection length ${L_{||}}$ itself is proportional to ${q_c}R$. Hence, $l_{||}^*$ does not change when varying ${q_c}R$, and we can assume it to be constant. Using parameters GA and GZ, the above equation can be written as

$$\begin{align}\lambda _q^{HD} \propto P_{SOL}^{1/8}{{q_c^{9/8}} \over {{B^{7/8}}}}G_A^{7/16}G_Z^{1/8}.\end{align} \tag{ A6 }$$

Substituting ${P_{SOL}}$ from equation (A1) into the above equation and replacing the GA dependence with the more relevant ${G_{AA}}$ dependence equation (4), as discussed in section 2, gives

$$\begin{align}\lambda _q^{HD} \propto {A^{0.31}}{{q_c^{9/8}} \over {{B^{0.78}}}} \times n_{e,sep}^{0.09}G_{AA}^{7/16}G_Z^{1/8}.\end{align} \tag{ A7 }$$

Here, the multiplication sign '×' separates the most important dependencies from the weaker dependencies.

In setting up new cases it is impossible to match the $\lambda _q^{HD}$ scaling in the code case, and potentially large correction factors ${f_{HD}}$ are inevitable. They may, however, be partly reduced by varying perpendicular transport coefficients, aiming at approaching the scaling

$$\begin{equation}{\lambda _{q,setup}} \propto {A^{0.31}}{{q_c^{9/8}} \over {{B^{0.78}}}}.\end{equation} \tag{ A8 }$$

To satisfy equation (A8), electron perpendicular heat conduction ${\chi _e}$ has to be varied in each code case. According to [23], equations (5.62) and (5.66), the relation between ${\chi _{e \bot }}$ and ${\lambda _q}$ obeys the scaling ${\lambda _q} \propto \chi _{_{e \bot }}^{1/2}$. Apart from this dependence, there is also the square root dependence of ${\lambda _q}$ on the product of upstream T_e and n_e , which is more difficult to include before knowing the output of the code run. In the present EDGE2D-EIRENE modelling, for the purpose of setting up individual cases, only the dependence ${\lambda _q} \propto \chi _{_{e \bot }}^{1/2}$ was used, giving for the ${\chi _{e \bot }}$ setup values

$$\begin{equation}{\chi _{e \bot ,setup}} \propto {A^{0.62}}{{q_c^{2.3}} \over {{B^{1.56}}}}.\end{equation} \tag{ A9 }$$