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

Predictions of the Huber–Chankin (HC) scaling for the upstream impurity fraction were verified in a series of EDGE2D-EIRENE (‘code’) runs for highly radiating plasmas with nitrogen injection. The main quantity extracted from the code was poloidally averaged, from X-point to X-point, separatrix impurity fraction cZ in the main scrape-off layer (SOL). Variation of the main working gas (H, D and T) revealed a qualitative agreement between the model and code results owing to the very large isotope difference in the predicted cZ values caused primarily by the inverse isotope mass dependence of the H-mode power threshold assumed in the HC model and implemented in the code. At the same time, the variation of the toroidal field and safety factor in deuterium cases yielded no correlation between the model predictions and code results. The code showed much higher local impurity fractions (fZ ) in the divertor compared to the main SOL, as well as large case-to-case variations in the divertor to the main SOL ratio of impurity fractions. The analysis of code results has wide-ranging consequences not only for the HC model, but also for other similar 1D models which use simple geometry ignoring strong neutral recycling in the divertor/ Different topology makes plasma parameters in the divertor and main SOL very different, resulting in different impurity charge state composition. Missing mechanisms in 1D codes (e.g. friction and thermo-forces exerted on impurity ions by main working gas ions) lead to impurity density redistribution. Neglecting all above factors, 1D models assume a constant impurity fraction along field lines.


Introduction
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 λ q is related to the experimentally established Eich scaling [9].In a simple 1D SOL (and divertor) model, λ 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 Hmode 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 λ q [11] instead of Eich's model, which does not have isotope mass dependence.The HD model for λ 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][5][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][14][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 thermoforces 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.

HC model and strategy of its verification by EDGE2D-EIRENE
The final scaling for c Z derived in [6] (equation ( 37)) is replicated here without a numerical coefficient and with some rearrangement of terms: 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, , 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 ne 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): In the above scaling, in its second proportionality, we replace ne with n e.sep and drop the f sep = n e,sep /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: where Ā is the average ion mass, including all ionization stages of impurities, 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 ε = 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 = λ q /λ HD q , where λ q is the power decay length equal to 2/7 λ Te , and λ HD q is given by equation ( 6) of [11].Since it is difficult to match λ q in the code calculations to λ HD q , the equations in [6] were re-derived to track the dependence on this additional parameter, yielding the dependence f −1.43   HD 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, Ā ≈ A and Z ≈ 1 (for hydrogen isotopes as the main gas), and the dependence of c Z on G A scales as G A ∝ A against A, resulting in the scaling c Z ∝ A −1.38 becoming much stronger: 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 ∝ 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 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: where 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 λ HD q as the power decay length, the exact match λ q = λ HD q 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 fcoefficients.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.

Setup of EDGE2D-EIRENE cases
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 2 ) 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 2 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.03MW 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.
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.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 λ HD q , which also has dependence on A.
The decay length λ 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][5][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 λ q = 2/7 × λ Te , where λ Te is the radial decay length of the electron temperature.In order to be able to change λ q in the code, the electron perpendicular heat conduction χ e has to be varied according to the scaling equation (A9) derived in appendix.The factor of the variation in χ 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 χ i ⊥ , but this procedure was not extended to the particle diffusion coefficient D ⊥ since, according to models [4][5][6], the n e decay length has no impact on c Z scaling, and only λ Te is taken into account when defining the SOL (power) width.
Generally, using variable χ ⊥ coefficients is not absolutely necessary, since any deviation of λ q from λ HD q extracted from the code output can be included in the correction factor f HD .
Still, varying χ ⊥ makes the code case parameters and profiles closer to those of the used model, reducing the correction factor f HD .

c Z from EDGE2D-EIRENE cases and predictions of the HC model
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 λ 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 ∝ n 0.72 e,sep 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 ∝ 1/n 0.72 e,sep , 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 ∝ 1/f 0.72 ne , giving the following total dependence on this parameter (and remembering that c Z scales with f 1.38  LH ): 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 95 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.
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.
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.
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.

Power balance in EDGE2D-EIRENE cases and related parameters
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.
f * LH is the normalized (on its value in the D case) coefficient fLH , 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, fLH is proportional to the ratio of P rad,N over P in : The dependence of c Z on fLH (or 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 where f corr is given by equation ( 6) as before.

Consistency between the full c Z scaling from the HC model and EDGE2D-EIRENE output
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.λ Te is calculated assuming exponential T e decay at the OMP position between poloidal rings s01 and s05 separated by 4.4 mm.λ 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 Ā, 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 Ā 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 Ā value highest in this case.Also, the average ion charge 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 λ 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 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.

Causes of discrepancies between the model and code c Z
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 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.
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.
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 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).
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-tocase 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 ∼λ q from the separatrix, better representing the middle of the 1D 'grid', since its width is determined by λ q in the HC model.The average λ 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%.

Summary
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: 'Ionizationdriven 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.

Figure 1 .
Figure 1.OMP profiles of Te, T i and ne for the reference D case.Vertical dash-dotted line indicates the separatrix position.

Figure 2 .
Figure 2. Target profiles of power density, Te and ne for the reference D case.Vertical dash-dotted lines indicate strike point positions.

Figure 3 .
Figure 3. Impurity fractions vs cell number for cases D, HB and LB.Vertical dash-dotted lines indicate entrances to divertor.

Figure 4 .
Figure 4. Te and ne profiles along field lines for cases D, HB and LB.Vertical dash-dotted lines indicate entrances to divertor.

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

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

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

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 fcorr, on the comparison between c Z values in the HC model and the EDGE2D-EIRENE code.

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