Validation of EDGE2D-EIRENE and DIVIMP for W SOL transport in JET

Tungsten sputtering rates and density profiles predicted using the edge plasma codes EDGE2D-EIRENE and DIVIMP are found to agree within a factor of 4 with measurements of neutral and singly-ionized W spectral line emission in the JET low-field side (LFS) divertor, and within a factor of 2 with SXR, VUV, and bolometric calculations of the W density in the main plasma. The edge plasma W predictions are extended to the core plasma using JINTRAC integrated core-edge modelling. Prompt redeposition of W is identified as the primary reason for the discrepancy between predicted and measured W emission in the divertor. The studied plasmas include attached divertor conditions in L-mode and type-I ELMy H-mode plasmas typical for JET. To more accurately reproduce the spectroscopically inferred W sputtering rates in EDGE2D-EIRENE, imposing the experimentally observed Be concentration of order 0.5% in the divertor is necessary. However, the W density in the main plasma is predicted to be insensitive to whether or not W is sputtered by Be at the divertor targets. Instead, the majority of the predicted core W originated in L-mode from sputtering due to fast D charge-exchange atoms at the W-coated tiles above the LFS divertor, and in H-mode due to D and W ions at the targets during ELMs.


Introduction
Tungsten (W) is the material of divertor plasma-facing components in ITER and the JET ITER-like wall [1], and a strong candidate for use in future fusion reactors.The erosion and transport of tungsten from the divertor plates into the core plasma is detrimental to the fusion performance due to the high radiation factor of W ions in the keV-MeV range [2].Therefore, the design of viable fusion reactors greatly benefits from the ability to predict the W erosion rate, the W screening by the scrape-off layer (SOL) and the main plasma W density as a function of time.
This study uses recent calculations [3] based on JET W diagnostics to quantitatively assess the ability of the 2D edge plasma codes EDGE2D-EIRENE [4,5] and DIVIMP [6,7] to predict the spatially resolved W sputtering and the main plasma W density in JET.The electron density and temperature profiles at the outer mid-plane as well as the beryllium densities in the divertor are adjusted to reproduce data from existing JET low-confinement (L-mode) and high-confinement mode (H-mode) experiments.The transport of W across closed flux surfaces is modelled by coupling the JETTO and SANCO core transport * Corresponding author.E-mail address: henri.kumpulainen@aalto.fi(H.A. Kumpulainen). 1 See the author list of E. Joffrin et al. 2019 Nucl.Fusion 59 112021.codes [8] with EDGE2D-EIRENE.All information from the W diagnostics is reserved solely for validation and not used as input in the W transport simulations.
The numerical tools used in this work have already been applied to several predictive studies on W transport in the edge and core plasma of JET [9][10][11][12] and ITER [12][13][14].The validation of these tools for W transport serves to investigate the accuracy and the regime of validity expected of such predictions, and to reveal potential areas of improvement.

Methods
One L-mode and one H-mode discharge (Table 1) with steady magnetic field, plasma current, heating power and electron density were selected for this study.Discharges with a horizontal low-field side strike point configuration (Fig. 1) were chosen due to the availability of high fidelity Langmuir probe and spectroscopic data for validation.While this configuration is not typical of the highest-performance JET baseline and hybrid scenarios, the selected discharges are similar to several https://doi.org/10.1016/j.nme.2020.100866Received 31 July 2020; Received in revised form 21 November 2020; Accepted 29 November 2020  other JET experiments.The focus of this study is on attached divertor conditions, due to the prevalence of W erosion in this regime and due to the limited reliability of the Langmuir probe data in detached plasmas.
In addition to the divertor conditions, also the upstream electron density and temperature were fitted to measurements by high resolution Thomson scattering [15], lithium beam [16] and reciprocating Langmuir probe [17] diagnostics when available (Fig. 2, 3).The fitting was done by adjusting the imposed values of ad-hoc cross-field transport coefficients.

EDGE2D-EIRENE simulation setup
The multi-fluid edge plasma/kinetic neutral code EDGE2D-EIRENE is used in this work to model the SOL plasma conditions, including intrinsic beryllium and tungsten impurities.The source of W is the predicted sputtering of divertor W surfaces due to D, Be and W. The Be concentration is based on the predicted sputtering of beryllium surfaces in the JET main chamber only; Be deposition and re-erosion is not self-consistently simulated.The sputtering model assumes Maxwellian ions, smooth surfaces with no roughness effects, and a perpendicular incidence angle.The 74 W ion charge states were bundled [18] to 6 fluid species for computational stability.
The applied heating power and fuelling rate were based on the experiments, taking into account the power losses due to radiation in the core plasma [19].For calculating the pumping rate, an effective albedo of 0.8 was specified for neutral particles at the LFS pumping plenum leading to the subdivertor.The HFS pumping plenum was assumed ineffective.These assumed albedo values are justified by earlier analysis of EDGE2D-EIRENE simulations with the subdivertor module included [20].
The anomalous cross-field transport coefficients for electrons and deuterium ions in the main SOL were selected to reproduce the experimental electron density and temperature profiles measured by highresolution Thomson scattering (HRTS) along the LFS mid-plane (Fig. 2).In the divertor SOL, the transport coefficients were set to 1 m 2 /s.For impurities, a constant  ⟂ of 1 m 2 /s was used in both the divertor and the main chamber.The pinch velocity was assumed 0 for all species.
In the H-mode scenario, ELMs are simulated with a frequency of 28 Hz and parameters based on earlier work [21].During the ELM phases, the particle and thermal diffusivities are multiplied by a factor of 10 and 30 respectively within an ELM-affected region on the LFS, extending radially from 0 to 7 cm inside the separatrix.Time-dependent kinetic correction factors [22] for the heat flux limiters were included, however corrections for the sheath heat transmission coefficients are neglected due to their adverse impact on the simulation time step when W is included.The duration of each ELM is specified as 1 ms.The recycling coefficient was assumed constant, neglecting any dynamic changes in the outgassing of the PFCs.Although this description of ELMs is very simplified and only applicable when the parameters are fitted to experiment, the goal here is to simply obtain background plasmas for W transport rather than attempting to predict the detailed ELM characteristics.The ELM parameters used in this study have been validated against statistical analysis of 53 similar JET type-I ELMy Hmode discharges with 12 MW of NBI heating power [21].Comparison of simulations with and without ELMs revealed that the increase in W erosion due to ELMs constitutes 80% of the predicted W in the main plasma, highlighting the importance of constraining the ELM parameters with accurate measurement data.This is consistent with earlier studies [23].
Test simulations were executed for L-mode and steady-state inter-ELM H-mode with cross-field drifts included.However, the drifts did not significantly improve the agreement for the primary target conditions relevant for sputtering.Furthermore, as drifts would not be compatible with the dynamic ELM simulations, due to numerical reasons, they were not included in this study.Generally, in other scenarios beyond this work, the drifts may be instrumental to reproducing the observed divertor conditions.The EDGE2D-EIRENE version used was v191219.

DIVIMP simulation setup
One significant advantage of the Monte Carlo code DIVIMP compared to the multi-fluid code EDGE2D-EIRENE is that all of the ionized states of tungsten can be individually tracked without adversely affecting the execution time and the stability of the simulation.DIVIMP has been observed to predict consistently higher W densities, by approximately 50%, than EDGE2D-EIRENE in JET L-mode and H-mode scenarios, primarily due to the bundling of W charge states in EDGE2D-EIRENE [24].Furthermore, EDGE2D-EIRENE simulations with both tungsten and cross-field drifts included tend to be unstable whereas in DIVIMP the inclusion of drifts has a negligible computational cost.These factors and the short execution time of DIVIMP make it an attractive option to use EDGE2D-EIRENE solely for the background plasma while replacing the computationally more demanding multifluid treatment of W with DIVIMP predictions, especially when drifts and/or ELMs are involved.
In this work, the background plasma conditions and the 2D neutral W ionization profiles are imported to DIVIMP from the EDGE2D-EIRENE simulations.Singly-ionized W test particles are launched at initial positions sampled from the W ionization profile.Non-prompt redeposition and thermal forces are enabled in both codes.The code version, boundary conditions and other options used are the same as in earlier comparisons between the two codes [24].

Integrated core-edge modelling setup
The experiment-based calculations [3] of the main plasma W density in JET are validated against soft X-ray diagnostics with error estimates in the plasma region with normalized minor radius < 0.6.Beyond this radius, the W densities are extrapolations based on bolometry.The SOL transport codes used in this work are not well suited to predict W transport in the core, hence the gap between the edge and core W densities is bridged by coupling EDGE2D-EIRENE to the 1.5D core transport codes JETTO and SANCO, the full system of codes collectively referred to as JINTRAC [8].The neutral deuterium atoms in the JETTO domain are described by a model called FRANTIC [25] and neutral beam heating is modelled using the PENCIL code [26].
The coupling interface between the EDGE2D-EIRENE and JETTO domains in the L-mode simulation is placed at a flux surface located 6 cm inside the separatrix at the LFS mid-plane.Coupling at the separatrix, while better justified by the transport models, excludes W erosion by D atoms from the confined plasma, because the kinetic distribution of neutral particles is not transferred across the interface.Including the 6 cm wide confined plasma region in EDGE2D-EIRENE is sufficient to reproduce the majority of the W sputtering by atoms observed in the standalone EDGE2D-EIRENE case.In the H-mode simulation, the interface is placed at the separatrix, because W erosion by D atoms has a low impact on the main plasma W density compared to the impact of the ELMs.
The neoclassical and anomalous cross-field transport is predicted using the standard NCLASS [27] and Bohm/gyro-Bohm transport models [28] respectively.The transport coefficients predicted by the core transport codes are also used as a core boundary condition in EDGE2D-EIRENE.In the H-mode simulation, a 2 cm wide edge transport barrier is set up in JETTO with local anomalous transport multipliers of 0.1, 0.03 and 0.025 for the electron thermal, ion thermal, and particle (green markers) and EDGE2D-EIRENE predictions with only intrinsic Be directly from the main chamber (black lines) and with 10 20 Be atoms/second injected at the targets to emulate Be re-erosion (orange lines).The JET divertor geometry is shown in red and the spectrometer lines-of-sight in blue.The predicted line emission was calculated using the pyproc post-processor [29].The excitation rates are from ADAS [30] with year identifier 89 for Be and 42 for W. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)diffusivities respectively.ELMs are modelled by JETTO, using the same transport multipliers as in the standalone EDGE2D-EIRENE case but with a Gaussian-shaped 5 cm wide ELM region centred at normalized radius 0.97.This approach was tested to produce similar heat and particle loads at the targets as the EDGE2D-EIRENE approach.No JETTO transport barrier or ELMs are applied to the integrated L-mode scenario.

Results
If EDGE2D-EIRENE is used to predict the Be flux to the LFS target by considering erosion and transport from Be components in the main chamber only, the resulting Be II line emission is lower than measured by a factor of 200 (Fig. 4(a), black line).The primary reason for this result is that the plasma-surface interaction model of EDGE2D-EIRENE neglects the re-erosion of Be from W surfaces and all of the deposited Be vanishes from the simulation.Fast reflection of Be is accounted for with coefficients based on TRIM data.Earlier studies have established that Be ions are the dominant cause of W erosion in typical attached JET Lmode scenarios and that a 0.5% concentration of Be 2+ at the LFS target is consistent with the observed W sputtering rate [31].Simulations with an unrealistically low Be concentration predict W erosion to be more localized at the strike point than observed by spectroscopy (Fig. 4(b) and (c), black line).This is due to thermal D ions only exceeding the W sputtering threshold near the temperature peak at the separatrix.The ion impact energy is affected by sheath acceleration and ion temperature.
The lack of a Be re-erosion model in EDGE2D-EIRENE can be compensated by injecting additional neutral Be at the targets to achieve the expected 0.5% concentration.In this case, both the predicted and observed Be II line emission profiles show a peak of similar shape, with the predicted Be II intensity 25% higher than the measurement (Fig. 4(a), orange line).The predicted W I emission is a factor of 3 lower than observed (Fig. 4(b), orange line), but the predicted peak W II emission is overpredicted by a factor of 4 (Fig. 4(c), orange line).The width of the predicted W I and W II emission peaks matches the measurements only when a realistic Be concentration is imposed.
Accurately reproducing both the W I and W II lines would require accounting for the prompt redeposition of W, for example by tracking the first few gyro-orbits of W ionized near surfaces.Prompt redeposition is expected for >90% of the W eroded under the studied plasma conditions [31], which is consistent with the discrepancy between the predicted and measured W I and W II emission in Fig. 4(b) and Fig. 4(c).One of the reasons why EDGE2D-EIRENE underpredicts the W I emission is the assumption of a perpendicular incidence angle for all ion species, which causes EIRENE to underestimate the sputtering yields.Surface roughness is another potential source of uncertainty.
The four rightmost spectrometer lines-of-sight intersecting with the outer vertical divertor indicate a second peak in the W I and W II emission, approximately one order of magnitude weaker than the peak emission at the strike point, despite the several orders of magnitude lower thermal load and incident ion flux in the far SOL.The EDGE2D-EIRENE simulations predict that the W erosion in this region is dominated by charge-exchange D atoms from the confined plasma and that this region is the largest contributor to the W influx across the separatrix in L-mode and inter-ELM scenarios.In contrast, the W erosion rate by Be and D ions at the strike points is predicted to have a negligible impact on the main plasma W density (Fig. 5, black and orange lines) due to efficient W screening by the SOL, except during ELMs when the W screening is weaker (Fig. 7, purple and orange lines).
The flux-surface averaged main plasma W density predicted by integrated modelling agrees with the experiment-based calculation [3] within 30% throughout the confined plasma.(Fig. 5).Such a close agreement can be attributed to coincidence, because even a 10%-15% uncertainty in the upstream SOL electron density would be enough to cause a 50% deviation in the predicted W density.Despite the sensitivity to the input parameters, the EDGE2D-EIRENE and DIVIMP predictions lie approximately within a factor of 2 from the extrapolated experimental result, demonstrating that they provide a viable boundary condition for predicting the main plasma W density.The 60% difference between the DIVIMP and EDGE2D-EIRENE predictions is mostly due to the bundling of W charge states to 6 fluid species in EDGE2D-EIRENE [24].
The toroidal rotation of the plasma causes a centrifugal effect which creates an asymmetry of up to one order of magnitude between the HFS and LFS main plasma W densities (Fig. 6).When a toroidal velocity of 50 km/s, consistent with JET measurements, is applied as a boundary condition at the core boundary in EDGE2D-EIRENE, the asymmetry is reproduced to within a factor of 2 on the closed flux surfaces (−1 < normalized radius < 1).DIVIMP does not include the centrifugal effect, thus predicting an opposite asymmetry driven by parallel-B SOL transport.The impact of rotation on the W density is only significant on the closed flux surfaces, due to parallel-B temperature gradients and inter-species friction playing a more dominant role in the SOL.Fig. 5. Flux-surface averaged W density in L-mode JPN 81472 at 9 s, calculated from diagnostics (red dashed line with shaded 30% confidence intervals), compared to code predictions.Orange solid and dotted lines: EDGE2D-EIRENE and DIVIMP respectively with 0.5% Be concentration at targets.Blue solid line: SANCO output from integrated modelling.Black solid line: EDGE2D-EIRENE without Be injection at targets.Green solid line: EDGE2D-EIRENE assuming zero toroidal plasma rotation.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)Fig. 6.W density along the mid-plane in L-mode JPN 81472 at 9 s.Negative radii correspond to the HFS and positive to the LFS.Linear interpolation is used to connect the EDGE2D-EIRENE and DIVIMP predictions across the core region.See Fig. 5 for legend.
The experimental calculation outwards of normalized radius 0.6 is extrapolated and expected to be valid until the pedestal top, but not in the SOL.The W density profile predicted by the core transport codes is one-dimensional, hence poloidal W asymmetries are not reproduced by the integrated modelling scenarios studied in this work.
The ELMy H-mode integrated modelling scenario overestimates the experimental W density by a factor of 2 near the magnetic axis and by a factor of 3 to 4 at normalized radii from 0.2 to 0.6 (Fig. 7).
Since it predicts a much larger W density gradient across the pedestal region than the edge codes and the measurements, the overestimation appears to be a feature of the core transport models in conjunction with the imposed ELMs and the edge transport barrier, rather than overestimated W influx by EDGE2D-EIRENE.More realistic integrated predictions could presumably be obtained by further constraining the JETTO pedestal parameters against experimental data.The standalone EDGE2D-EIRENE and DIVIMP W density predictions are intended to be valid for the SOL, but of limited use far inside the separatrix due to the purely diffusive radial transport model.

Conclusions
When the background plasma conditions are fitted to measured JET L-mode and type-I ELMy H-mode plasmas in attached outer divertor conditions, the plasma codes and DIVIMP demonstrated ability to predict the neutral and singly-ionized W emission in LFS divertor within a factor of 4 and the upstream W density to within a factor of 2 of the experimental values.One of the probable reasons to the underestimated neutral W emission is that EDGE2D-EIRENE assumes a perpendicular incidence angle for all ion species.On the other hand, the singly-ionized W emission is overestimated due to not properly accounting for W prompt redeposition.The validity of the main plasma W predictions is supported by integrated core-edge modelling, using JETTO and SANCO for core transport.Keeping in mind that the codes used in this work involve simplifications which are not necessarily accurate under all circumstances, further validation is necessary when extrapolating these results to different plasma conditions.
The uncertainties resulting from the simulation parameters in predicting the W density are greater than the observed discrepancy between the predictions and the experiments.The predicted main plasma W density is particularly sensitive to the upstream electron density, the ELM transport multipliers and the width of the ELM-affected region.For predicting the W density in future discharges, it is thus very beneficial to constrain these parameters as accurately as possible using the available data from existing discharges as done in [21].
Accounting for Be re-erosion at the targets greatly affects the predicted W erosion rate, but has little effect on the main plasma W density due to effective divertor screening.Including toroidal rotation of the plasma in the simulations was shown to reproduce the experimentally observed HFS-LFS W asymmetry due to the centrifugal effect, without affecting the W density in the SOL.
Diagnostic uncertainties in the upstream electron density and temperature profiles, of the order of 15%, are alleviated by simultaneously validating the target conditions which are highly sensitive to the upstream values.Uncertainties in the comparison of synthetic and measured W I and W II emission result from model simplifications as well as the uncertainties in the plasma conditions and sputtering and emission data, cumulatively at least a factor of 2. While the values of the crossfield transport coefficients are uncertain, their impact on the predicted W density is considerably weaker than a linear dependence.For studies assessing the effects of certain model assumptions, such as alternative sputtering models, surface mixing of Be and W, or prompt redeposition, a different choice of simulation tools may be more appropriate.

Table 2
Catalogue identifiers of the EDGE2D-EIRENE and integrated core+edge JINTRAC simulations stored on the JET Heimdall cluster (July 2020).The DIVIMP simulations with corresponding identifiers are stored on the IPP Garching TOK cluster and on the Aalto University Triton cluster.

Fig. 1 .
Fig. 1.Divertor target configuration and discretized computational domain of JPN 81472 at 49 s.The SOL is shown in blue, the main plasma region in purple, and the private flux region in green.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 2 .
Fig. 2. Ion and electron cross-field heat and particle diffusion coefficients (a, d) used in EDGE2D to reproduce the electron density and temperature profiles measured by high-resolution Thomson scattering, Li-beam, and reciprocating probe along the LFS mid-plane (b, c, e, f) in the L-mode (a-c) and inter-ELM H-mode (d-f) scenarios.

Fig. 3 .
Fig. 3. Electron temperature (a, c) and density (b, d) profiles predicted by EDGE2D-EIRENE (solid lines) compared to Langmuir probe (LP) measurements (markers) along the LFS target in the L-mode (a, b) and inter-ELM H-mode (c, d) scenarios.The predicted target electron temperatures approximately match the measurements, whereas the predicted electron densities exceed the LP measurements by approximately 30% in L-mode and by a factor of 2-3 in inter-ELM H-mode.The LP data are limited to a maximum of 100 eV and thus are not shown for the intra-ELM conditions.

Fig. 4 .
Fig. 4.Comparison of (a) Be II, (b) W I, and (c) W II spectral line emission from JPN 81472 at 9.5-10.0s measured by the high-resolution Czerney-Turner spectrometer KT3 (green markers) and EDGE2D-EIRENE predictions with only intrinsic Be directly from the main chamber (black lines) and with 10 20 Be atoms/second injected at the targets to emulate Be re-erosion (orange lines).The JET divertor geometry is shown in red and the spectrometer lines-of-sight in blue.The predicted line emission was calculated using the pyproc post-processor[29].The excitation rates are from ADAS[30] with year identifier 89 for Be and 42 for W. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 7 .
Fig. 7. Flux-surface averaged W density in H-mode JPN 82486 at 14 s, calculated from diagnostics (red dashed line with shaded 30% confidence intervals), compared to code predictions.Orange solid line: EDGE2D-EIRENE with ELMs.Blue solid line: SANCO output from integrated modelling with ELMs.Purple solid and dotted line: EDGE2D-EIRENE and DIVIMP respectively without ELMs.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) H.A. Kumpulainen: Formal analysis, Investigation, Visualization, Writing -original draft, Data curation.M. Groth: Supervision, Conceptualization, Writing -review & editing, Project administration.G. Corrigan: Methodology, Software.D. Harting: Methodology, Software.F. Koechl: Methodology, Software.A.E. Jaervinen: Software, Conceptualization, Writing -review & editing.B. Lomanowski: Methodology, Software.A.G. Meigs: Methodology, Software.M. Sertoli: Methodology, Formal analysis.

Table 1
The main parameters of the studied JET discharges JPN 81472 at 9 s in L-mode and JPN 82486 at 14 s in H-mode.