Full-radius integrated modelling of ASDEX Upgrade L-modes including impurity transport and radiation

An integrated framework that demonstrates multi-species, multi-channel modelling capabilities for the prediction of impurity density profiles and their feedback on the main plasma through radiative cooling and fuel dilution is presented. It combines all presently known theoretical elements in the local description of quasilinear turbulent and neoclassical impurity transport, using the models TGLF-SAT2 and FACIT. These are coupled to the STRAHL code for impurity sources and radiation inside the ASTRA transport solver. The workflow is shown to reproduce experimental results in full-radius L-mode modelling. In particular, a set of ASDEX Upgrade L-modes with differing heating power mixtures and plasma currents are simulated, including boron (B) and tungsten (W) as intrinsic impurities. The increase of predicted confinement with higher current and the reduction of core W peaking with higher central wave heating are demonstrated. Furthermore, a highly radiative L-mode scenario featuring an X-point radiator (XPR) with two intrinsic (B, W) and one seeded argon (Ar) species is simulated, and its measured radiated power and high confinement are recovered by the modelling. The stabilizing effect of impurities on turbulence is analysed and a simple model for the peripheral X-point radiation is introduced. A preliminary full-radius simulation of an H-mode phase of this same discharge, leveraging recent work on the role of the E×B shearing at the edge, shows promising results.


Introduction
The ability to predict the behaviour of impurities has become critical as magnetic fusion devices approach reactor-relevant a See Zohm et al 2024 (https://doi.org/10.1088/1741-4326/ad249d)for the ASDEX Upgrade Team.operation.High-Z metallic environments are planned for next generation tokamaks [1][2][3] and material limitations on incident excess heat fluxes will likely require high impurity seeding for safe power exhaust [4], with high core radiative power fractions (f rad ≡ P rad /P heat > 50%) envisioned for reactor scenarios [5][6][7][8].At the same time, core contamination with heavy impurities can lead to unacceptable radiative losses of plasma energy, while light impurities can dilute the fuel to insufficient levels for economical fusion power generation [9].
Full-radius integrated modelling capabilities for accurate predictions of confinement have recently been demonstrated in tokamaks, both in H-mode (where the complex pedestal region requires a description of not only transport but also magnetohydrodynamic (MHD) stability) [10][11][12][13] and L-mode (where an accurate description of turbulent transport up to the separatrix is necessary) [14,15].In particular for the latter, the TGLF-SAT2 quasilinear transport model [16] was shown to reproduce ASDEX Upgrade (AUG) experimental results in full-radius simulations [14].The use of theory-based transport models allows for a more confident extrapolation to machine parameters which are unexplored by scaling law databases, and 1.5D integrated modelling provides not only global quantities but entire profiles of the plasma parameters.So far, however, these full-radius modelling frameworks do not selfconsistently describe impurity transport.
Recent developments in integrated modelling have been dedicated to the description of impurity transport [17][18][19][20][21][22][23][24], including the evolution of the background plasma with radiative losses by impurities [25][26][27].However, some limitations are present either because the boundary condition is set well inside the confined plasma, or because empirical impurity transport coefficients are used at the edge.To our knowledge, physics-based quasilinear models have not been applied in a full-radius description of impurity transport so far.Likewise, the effects of high impurity contents in strongly radiative scenarios deserve further consideration within integrated modelling efforts.
In this work we present an integrated modelling framework that self-consistently describes impurity transport and the feedback effects of impurities on the main plasma via radiation and fuel dilution.The framework is described in detail in section 2. It is then applied to the modelling of well-diagnosed AUG experiments in sections 3 and 4, where full-radius simulations of unseeded L-modes and a strongly seeded high confinement L-mode are presented, respectively.Section 5 provides a summary and outlook.

Integrated modelling framework with impurity transport
The models for transport, heat and particle sources, and equilibrium are integrated in ASTRA [28,29].The approaches for the main plasma (electrons and fuel ions) and impurities are discussed in the following subsections.The workflow is summarized in figure 1.

Main plasma
The modelling of the main plasma is mostly based on that of [14], and it is only briefly summarized here.
Turbulent transport coefficients for the evolution of T e , T i , n e are calculated with TGLF-SAT2, whereas NCLASS [30] is used for the neoclassical transport.We use 161 radial points in ASTRA, calling TGLF-SAT2 every second point.Neutral beam injection (NBI) and electron cyclotron resonance heating (ECRH) sources are calculated with RABBIT [31] and TORBEAM [32] respectively.The particle source is set by the density of neutrals entering the separatrix, with a feedback loop on the volume-averaged electron density ⟨n e ⟩ vol , in analogy to the active control of the line-averaged density n e in the experiments.⟨n e ⟩ vol is used in order to preserve the total particle content and, therefore, the heating power per particle.We assume the neutrals are Franck-Condon atoms with an incoming energy of E 0 = 2 eV and neglect charge exchange (CX) particle sources.This method to determine the particle source in the simulations has been tested in detail in [14].The main ion density is calculated from n e and the impurity densities, by imposing quasi-neutrality.The current density profile is also evolved, including the bootstrap current formulae from [33].The 2D magnetic equilibrium is calculated in ASTRA from the simulated pressure and current profiles using SPIDER [29,34], setting the separatrix shape obtained from CLISTE [35] reconstructions.
All profiles are simulated up to the separatrix.As boundary condition we use the 2-point model for T e,sep [36], with a heat flux width λ q from the heuristic drift model [37], and we set T i,sep = 1.5 T e,sep , n e,sep = 0.3⟨n e ⟩ vol like in [14].The radial electric field E r becomes particularly relevant at the edge of ion heated discharges at powers approaching the L-H transition, through its role on the E×B shearing [38].Our description is based on the work of [39], where the edge E r well is mimicked by solving the ion force balance until a radial position ρ min in the normalized poloidal flux ρ pol , then forcing E r to zero at the separatrix (thus E r has a local minimum at ρ min ), such that where α = 2 and we use ρ min = 0.985, informed by AUG Lmode measurements [40,41].This is a new element with respect to [14].Further discussions on this treatment of E r and a sensitivity study on the choice of ρ min are presented in section 3.3.

Impurities
The evolution of an impurity with flux surface average (FSA) density ⟨n z ⟩, diffusive and convective transport coefficients D z and v z , and sources S z is described by where V is the volume enclosed by the flux surface labelled by some radial coordinate r, V ′ = ∂V/∂r, and g r = ⟨(∇r) 2 ⟩ is a metric coefficient.The sources/sinks of a charge state are determined by atomic processes of ionization and recombination with the neighbouring ionization stages and CX reactions with the fuel [42].
The turbulent component of the impurity transport coefficients is calculated with TGLF-SAT2.It has been shown in [43] that there are non-negligible centrifugal effects on the turbulent transport of heavy impurities like tungsten (W), due to the increased impurity density on the low field side (LFS), where the turbulence balloons.These effects can be included in post-processing via the analytical formulae derived in [43], as done already in the QuaLiKiz model [44].We have now implemented this missing physics in the TGLF interface to ASTRA, under the simplifying assumptions that the poloidal functional form of the amplitude of the electrostatic potential fluctuations is given by | φ| 2 ∝ exp [−θ 2 /(2 θ rms ) 2 ], and that its mean width θ rms follows the GLF23 parametrization [45] where q is the safety factor and s = (r/q) dq/dr is the magnetic shear.Careful transformations from LFS to FSA coefficients are necessary, and are done following the appendix of [46].The turbulent coefficients are enhanced by a factor of ⟨⟨exp(−E z )⟩⟩/⟨exp(−E z )⟩ ∼ 1-3, ranging within typical values of the rotation.Here, the FSA is denoted by single brackets ⟨ • ⟩, and the mode envelope average ⟨⟨ • ⟩⟩ is defined as where J is the Jacobian of the coordinate system.If the background electrostatic potential Φ(r, θ) is perturbed only by charge separation due to toroidal rotation, and not by temperature anisotropies introduced e.g. by ICRH [47], then the normalized impurity energy is where R is the major radius coordinate, with R LFS and R 0 its values at the outer mid-plane and geometric axis respectively.The effective impurity Mach number is and the main ion Mach number is the ratio of toroidal rotation and main ion thermal speeds, Neoclassical transport is calculated with NCLASS for low-Z impurities and FACIT [48][49][50] for high-Z impurities.The latter describes the effects of rotation-induced poloidal asymmetries at all collisionalities, which are essential for heavy impurities given that they greatly enhance neoclassical transport [51][52][53].The applicability of NCLASS at high collisionalities is limited, as it overestimates the temperature screening of impurities [49], unlike FACIT and the more complete drift-kinetic code NEO [54][55][56].However, the collisionality windows of the applications in this work are still in the domain where these three models are consistent.In contrast, NCLASS describes multi-species collisions, while FACIT considers main ion-impurity collisions and neglects at present collisions between different impurities (an approach to overcome this limitation is proposed in section 4).Instead, the much stronger effects of rotation on heavy impurities are described by FACIT, which NCLASS cannot calculate.
The impurity code STRAHL [57] is used to calculate the ionization equilibrium of the impurity species from their sources and transport coefficients, as well as the radiated power profiles they produce.This allows us to consider the effect of impurities on the local power balance through radiative cooling of the electrons, in addition to the fuel dilution that arises by imposing quasi-neutrality.For the calculation of D z and v z in TGLF, FACIT and NCLASS, the input charge of the impurities at a given radial position is their average charge calculated by STRAHL, assuming that in the close vicinity to that location charge states far from the average charge are not very populated.STRAHL calculates the value of the impurity densities at the separatrix from the input edge sources and the simplified scrape-off layer (SOL) model described in [58].For this model we use typical values of the SOL Mach number [59], and we calculate the SOL decay lengths of the plasma profiles, λ, assuming a Spitzer-Härm conductivity relation λ Te = 7λ q /2 along with the experimental observation of λ ne ≈ 3λ Te /2 [60], and taking λ Ti ≈ λ Te , λ nz ≈ λ ne (for all impurity species).Collisional thermalization of the impurities, T z = T i , is assumed throughout the plasma radius.
The coupling of STRAHL and ASTRA has been updated and generalized to an arbitrary number of species.A scheme to correctly split the turbulent flux into diffusive and convective components has been implemented for TGLF, which does not output the individual transport coefficients but rather the total flux.This is essential for realistic ionization equilibrium calculations by STRAHL.A correct diffusion coefficient is particularly important due to the strong gradients of the individual charge state densities, even if the resulting average density profile has a mild gradient.The method to obtain such splitting consists in performing two TGLF runs, one with the real impurity density gradient and one setting the gradient to zero.In the latter the entire impurity flux is convective, Γ 2 /n = v, which can be replaced in the first computed flux, Here, the inverse gradient length scale of a profile X(r) is 1/L X = −∂ ln X/∂r.For this method we assume that the impurity gradients do not strongly affect the linear stability of the background turbulence.We choose this approach over another possibility of including an additional trace impurity with zero gradient for each simulated species, because the cost of the TGLF matrix inversion grows with the square of the number of species.With two background species (electrons and main ions), the computational times of the twocall and duplicate-trace methods grow as 2 × (2 + N imp ) 2 , (2 + 2 N imp ) 2 for N imp impurities.The two-call method is already faster for two impurities, the minimum we consider in our applications in sections 3 and 4. For parameters where the fluxgradient dependence is highly non-linear this method might not be suitable, but tests with our AUG parameters have not shown this to be the case.
We set a significantly lower time step for STRAHL than ASTRA (10 µs vs 10 ms), because STRAHL solves the impurity equations explicitly in time with ionization-recombination times much shorter than 10 ms.A lower time step is also required for numerical stability, particularly in cases with strong convection, given that the transport solver of STRAHL lacks the numerical scheme for stiff transport [61] of ASTRA.

Additional elements
An additional ion species, corresponding to the fast ion density n fi from RABBIT, is kept in ASTRA in the quasineutrality condition (but not included in TGLF as a kinetic species) to model fast ion dilution.Fast ions can stabilize turbulence [62,63] and can have a relevant influence on the transport of impurities [64].
Effective diffusivities in the inner core region (normalized toroidal flux ρ tor < 0.3) are introduced to account for sawtooth relaxations, with values between 0.1-1.0m 2 s −1 .This is done for the kinetic profiles in order to provide more stability to the full-radius simulations with TGLF, whereas the Kadomtsev reconnection model is applied to the safety factor and current density only.
Finally, note that toroidal rotation is an experimental input.Momentum transport is the main channel that we do not simulate at present, considering the difficulties in predicting the toroidal rotation profile accurately, particularly in conditions with low or zero externally applied torque.Therefore, selfconsistent momentum transport modelling falls outside the scope of this work.However, we reiterate that heavy impurity transport is sensitive to the rotation, both its magnitude (which affects both neoclassical and turbulent components, albeit more strongly the former) and its shear (which affects the turbulent roto-diffusion).

Full-radius simulations of AUG L-modes
We consider a set of six AUG L-mode phases in deuterium, which differ mostly in their NBI-ECRH heating mix and plasma current.The main parameters of each discharge are summarized in table 1.These experiments have no seeded impurities.We assume only W and boron (B) are present, entering the plasma via its interaction with the boronized tungsten walls.
In addition to diagnostics measuring the main plasma, namely the integrated data analysis suite (IDA) [65] for n e and T e , and charge exchange recombination spectroscopy (CXRS) for T i and v φ [66], CXRS measurements of n B [67] were available for all discharges except for the phase without NBI, as well as soft x-ray (SXR) data [68] that allow us to extract n w .Likewise, bolometric estimations of the total radiated power [69] and IDA-based estimations of Z eff (IDZ) [70] were available for comparisons to other diagnostics.

Impurity sources in the simulations
An important element of these simulations is how to set the impurity sources, which are the main input of STRAHL.For a given profile shape, determined by transport, these sources ultimately determine the content of each species by setting the value of its density at the separatrix.In STRAHL, this is determined by neutral puff rates for each impurity and the model for SOL losses.In the experiments neither B nor W are puffed, but rather enter the plasma due to wall erosion and sputtering.However, a consistent modelling of real sources is outside the scope of our work.To compare impurity transport predictions and measurements one is interested in the impurity density gradients, because R/L nz = −R v z /D z in steady-state, Table 1.AUG L-modes to be analysed in section 3. The time column is the interval (in s) over which profiles are averaged, the injected NBI and ECRH powers are given in MW and the plasma current Ip in MA. q 95 is the safety factor at 95% of the toroidal flux.AUG has an aspect ratio of 3.3 and R 0 = 1.65 m.The magnetic field and line-averaged density are similar for all cases: Bφ ≈ 2.5 T and ne ∈ {2.0−2.6}× 10 19 m −3 .The red and orange cells denote heating mix scans (at constant Ip), the blue and cyan cells denote Ip scans (at constant heating mix).whereas the total content of impurities is relevant in terms of radiative losses and fuel dilution.

Shot
If we estimate Z eff using only CXRS B densities, we find a consistent underestimation of the IDZ measurements.B is the most prominent light impurity in AUG plasmas due to routine wall boronizations, however several other species can be present in smaller concentrations to yield a non-negligible contribution to Z eff .In particular, nitrogen (commonly seeded for diagnostics and heat exhaust) and helium (from glow discharges performed for the boronizations) can remain in the vessel from previous experiments, with smaller traces of oxygen, carbon and fluorine.Simulating all these species is unaffordable.We opt for considering that only B is present, but in a higher content that includes the other light impurities, assuming that the transport of these species is not too different due to the similar charge and charge-to-mass ratio.Furthermore, we observe that the measured line-averaged effective charge Z eff is quite constant across these six phases, at Z eff ≈ 1.7.The approach to select the B source consists on a feedback loop on the B neutral puff rate such that the simulated Z eff reaches a target value of 1.7 for all simulations.The W source is set by feedback on the simulated total radiated power, such that the total radiation estimated by bolometry is matched.
Figures 2 and 3 show relevant profiles of two phases with different engineering parameters, AUG #39255 at 2.5-3.0 s (high current, mixed heating) and AUG #39323 at 5.0-6.0 s (low current, NBI heated), with particular interest on the contrast between their impurity and radiation profile shapes.Equivalent figures for the other four phases are not shown for brevity, but general properties of the full set of simulations are discussed in section 3.3.
The temperature profiles of the mixed heating discharge, shown in figure 2(a), are well reproduced, although they are somewhat overpredicted in the core.In figure 2(b), we see that the predicted n e has a weaker edge gradient and slightly higher core density.Recall that the experimental volume-averaged density is matched by a feedback on the neutral source, so if the transport does not predict a strong enough edge gradient the source increases, raising the entire n e profile.The simulated Z eff , in figure 2(c), matches the experimental value by design.The safety factor profile q of this high current case, shown also in figure 2(c), has an inversion radius at ρ tor ≈ 0.4, inside which the reconnection model keeps q = 1.The simulated n B , shown in figure 2(d), is peaked and has a lower edge gradient than the experiment, a similar gradient at mid-radius and a larger gradient in the core, where the experimental profile is flat.At the edge, the B 5+ density of the simulation is more consistent with the B 5+ measurement, since this charge state is less populated at the colder periphery.Following section 3.1, the B content that generates the target Z eff is 2.6 times larger than the CXRS measurements.The experimental W density is calculated from FSAs of 2D SXR tomographic reconstructions [71,72], assuming that only W generates the emissivity ϵ sxr (MW m −3 ) (Bremsstrahlung is negligible in these cases), so that where L sxr w (MW m 3 ) is the W SXR cooling factor.Note that the absolute calibration of the SXR signal is highly uncertain.We typically need a scaling factor that is estimated such that the volume integral of the radiated power density caused by ⟨n exp w ⟩ approximately matches the total radiation measured by bolometry.For instance, calculating the radiated power caused by a W density derived from the raw SXR data in AUG #39255 leads to approximately twice the radiation of the bolometry, so we multiply the SXR emissivity by 0.5, as shown in the legend of figure 2(e).Nonetheless, SXR data are extremely valuable in terms of the profile shape and gradients they provide and the comparisons to the simulated W transport they allow, particularly from axis to mid-radius.The SXR cameras have a lower detection limit at T e ≈ 1 keV [68], below which the data is unavailable, so we lack edge W profiles.
The measured core W profile in figure 2(e) is hollow, consistent with the presence of central ECRH.The simulated n w is flat and even slightly hollow in the core, due to the reduced neoclassical pinch caused by the flat n i close to the axis and a mild central toroidal rotation of v φ ≈ 78 km s −1 , M i ≈ 0.15.The simulated profile of the radiated power density shown in figure 2(f ) is therefore also quite hollow.Direct measurements of this quantity were not available, however from bolometry we have that the total radiated power from inside the separatrix (that is, the volume integral of the radiation density) is 1.2 ± 0.3 MW.The uncertainty of the bolometry is taken as the difference between the minimum and maximum of the signal within the time window under consideration.The contribution of B to the simulated total radiation is not negligible and it is localized at the very edge of the plasma, where the lower temperatures allow for line emissions of partially ionized B.
The main profiles of the low current, beam heated discharge are presented in figure 3.Both simulated temperature profiles are close to the experimental data, as shown in figure 3(a).Note the strong edge gradient in the experimental T i .The high core n e peaking is well reproduced, but its large measured edge  gradient is not, as shown in figure 3(b).The Z eff profile, shown in figure 3(c), becomes quite peaked deep in the core, where the W content is not negligible.The edge q profile is large as expected, and the inversion radius lays at ρ tor ≈ 0.2.Like n e , the measured n B of figure 3(d) has a strong edge gradient, but it is then flat until mid-radius after which it becomes slightly hollow and then quite peaked deep in the core.The simulated average B density does not capture these strong changes in the profile shape, being mildly peaked throughout the radius, however the simulated B 5+ agrees more with the steep edge gradient as before.In figure 3(e), the SXR data results in a very peaked experimental n w , whose profile shape is well reproduced by the simulation.The combination of a very peaked central n i and a high rotation (M i ≈ 0.36) gives rise to the strongly inward neoclassical convection that causes this peaking.The radiated power density in figure 3(f ) is therefore also very peaked, and the 0.5 ± 0.1 MW measured by bolometry are well matched.Note that the B contribution to the total radiation is even dominant in this case, given that W has accumulated in the core whereas B radiates in the edge, which has a larger share of the plasma volume.

Discussion
Some general properties of the set of simulations are presented in figure 4.An important initial test for the integrated modelling framework is to verify that the trend of increased confinement with higher current is reproduced.This is investigated in figure 4(a), where a clear current dependence of the stored thermal energy W th is observed and the experimental values are well matched by the simulations.It has been shown in [15] that including the E×B shearing rate γ E×B in ASTRA/TGLF-SAT2 simulations is an essential element to capture this dependence.This is calculated in our simulations with the Waltz-Miller formula [73] Note that the ITER physics basis (IPB) scaling of L-mode confinement is ∝ I 0.96 p [74].While the limited set of simulations hinders a proper statistical analysis, we can nonetheless fit the current dependence of the six simulations for comparison to scaling laws and previous modelling results.We recover a dependence of I 0.87 p that is consistent with the I 0.84 p of [14] and lower than the almost linear IPB scaling.Note that while the auxiliary power is approximately matched in these discharges, the total power is not, considering the different Ohmic heating.Accounting for this reduces the current scaling due to the small increase of stored energy produced by the Ohmic power at higher I p [14].
Figure 4(b) presents a correlation between the mid-radius electron and boron density gradients.The trends of the measurements and the simulations coincide.There is, however, a general overestimation of the n e peaking.Note the missing experimental B gradient in the purely electron heated discharge.There is no CXRS data with only ECRH, because the CX reactions occur between the beam ions and the impurities.This is unfortunate, considering it is the flattest n B found by the modelling.Standalone TGLF calculations have been performed to better analyse the difference in mid-radius turbulent transport of boron in two extreme cases at the same current (0.8 MA): a purely beam-heated phase with more peaked boron and electrons (orange square) and an ECRH-only phase with flatter boron and electrons (green square), corresponding to AUG #35475 at 2.8-3.5 s and 5.8-6.5 s respectively.The diffusive and multiple convective components are split by setting the corresponding B gradients to zero, following the general expression for the turbulent impurity flux in equation (20) of [75].The case with beam heating has a dominant inward pure pinch that is only partially compensated by smaller outward thermo-and rotodiffusive fluxes.In contrast, the ECRHheated case presents an outward pure pinch component that is slightly overcome by inward thermodiffusion (with negligible rotodiffusion), leaving a smaller inward convection.Finally, note that the CXRS profile at low current is hollow, but the modelling predicts no hollow profiles.The lack of predicted light impurity hollowness is a well documented missing piece of theoretical understanding [76][77][78][79], although two physical mechanisms have been recently proposed to explain it [80,81].
Figure 4(c) studies the effect of ECRH on the core W peaking. Applying central wave heating to control heavy impurity accumulation is an established experimental technique [82][83][84] with a solid theoretical understanding [85,86].In the set of simulated discharges we have a scan in ECRH power fraction at approximately constant total heating, which allows for comparisons between the SXR-measured and predicted ratio of W densities on and off axis.We find a strong reduction of the W peaking with increasing ECRH fraction, in agreement with the experiment.
We conclude this section by discussing the effect of the E r boundary condition, described in section 2.1.A sensitivity study on the choice of radial location for the minimum of the E r well is shown in figure 5, for the same discharge as in figure 3, where the measured T i data shows a strong edge gradient.Here, 'floating' means E r is allowed to go as negative at the separatrix as the simple ion force balance predicts.Bonanomi et al recently showed that TGLF-SAT2 can form pedestal-like structures [39], in situations where the simulation enters a feedback loop: γ E×B reduces the edge turbulence, which in turn raises the edge pressure gradient, deepening the E r well, thereby increasing γ E×B .This happens in figure 5 for the cyan and orange curves, which largely overshoot typical minimum E r values of −5 to −10 kV m −1 for L-modes at AUG [41].This leads to considerably stronger simulated edge pressure gradients than the experiment.However, decreased core gradients in these two cases partially compensate the increased edge pressure, such that the global confinement is 'only' overestimated by 20%.The magenta line corresponds to the simulation shown in figure 3, for which the E r well is instead somewhat underpredicted.The edge gradients are weaker but the core pressure is higher, and the confinement is only 5% higher than the experiment.This is the boundary condition chosen for all L-mode simulations in this paper.For reference, a simulation without E×B shearing in TGLF underpredicts both edge and core, leading to a 13% lower simulated confinement.
Given that the edge impurity convection is mostly neoclassical (so it is quite sensitive to the main plasma gradients) and the impurity turbulent diffusivity is also reduced by the E×B shearing, these elements are critical for the peripheral impurity transport.physics, such as the impact of self-generated Reynolds stresses on the different rotation terms entering in the radial force balance [87][88][89] or effects due to ion orbit losses [90][91][92], could in the future provide a self-consistent description of the E r well that is better reconciled with the positive SOL E r values set by parallel dynamics on the open field lines and sheath boundary conditions at the target [93].This would go beyond the artificial treatment presently used in equation (1), which mimics the experimental observations of the edge E r in L-mode conditions, but it is quite outside the scope of our work.

High confinement radiative L-mode with X-point radiation
The experiment presented in [94], AUG #37041 at 5.0-5.5 s, is a high-power (P NBI = 5 MW, P ECRH = 2 MW) radiative Lmode kept just below the L-H power threshold P LH by active feedback control on the power crossing the separatrix, P sep = P aux − P rad , using argon (Ar) seeding as actuator.It also features a region of strong radiation localized above the X-point.An interesting property of this discharge is that, after an initial H-mode phase, it transitions back to L-mode but retains high confinement (H 98 ≈ 0.95) without ELMs.This type of high performing radiative L-modes has been obtained in different machines [95][96][97][98][99] and is currently investigated for potential reactor scenarios [100].AUG #37041 is an appealing application due to the high impurity content and radiated power, the presence of multiple species (intrinsic B and W, seeded Ar) and the high confinement with no ELM activity.Furthermore, it is a well-diagnosed discharge, with CXRS measurements of Ar 16+ in addition to the diagnostics listed in section 3.2.On the other hand, it is also a challenging application since, despite being a stationary phase, the balance between heating and radiated powers is marginal, so a slow increase of the radiation can cause a radiative collapse (which in fact happens later on in the experiment).

Impurity sources in the simulation
The approach to determine the impurity sources in STRAHL for the simulation of this radiative L-mode differs from the one discussed in section 3.1.Here the seeded argon has dominant contributions to both the radiated power and the effective charge.The neutral Ar puff rate in STRAHL is increased proportionally to (P sep − P LH ), calculating P LH with the Martin scaling [101], in close analogy to the experiment.The W source also follows this proportionality, although with a feedback loop 10 times weaker, with the idea of mimicking a higher W sputtering by a higher number of Ar particles.In contrast, the B source is set as a constant trying to match the line average of the CXRS-measured B profile.

Modelling results against experimental data
Figure 6 presents relevant profiles of this experiment.T e , T i and n e are well reproduced (particularly the latter two, with a somewhat overestimated core T e ), as shown in figures 6(a) and (b).In figure 6(c), the high Z eff of the discharge, approaching 3 due to the high Ar seeding, is matched within uncertainties.The q profile has an inversion radius at ρ tor ≈ 0.3.The core n B is more peaked than the flat CXRS profile, as shown in figure 6(d).The peripheral transport reproduces the gradient at ρ tor > 0.7, but we lack CXRS data at ρ tor > 0.86 for complete edge comparisons.
In this experiment, extracting a W density from SXR measurements is not as simple as in section 3, because the high Ar content means that it dominates the SXR emissions.
For the W profile, two simulations are run with different neoclassical transport from FACIT: the first considers collisions with the main deuterium ions only, and the second has an additional flux component where FACIT is applied with argon as the main ion.The latter approach has not been systematically validated with respect to NEO (unlike the former, in [49,50]), but for the specific profiles in this simulation, standalone NEO calculations show quantitative agreement with the increase in the magnitude of the neoclassical diffusion and qualitative agreement with the decrease of the temperature screening associated with a higher effective collisionality when collisions with Ar are considered.Both resulting W profiles are shown in figure 6(e).The profiles are quite similar for most of the radius, being slightly peaked in 0.3 < ρ tor < 1. Inside ρ tor ≈ 0.3 both profiles peak considerably, but the one without Ar collisions then flattens.The increased core peaking in the case with Ar collisions is due the more inward neoclassical convection associated with the reduced temperature screening.All other profiles in figure 6 are not noticeably modified between the two simulations, so moving forward we consider only the simulation with collisions of W and D only (where FACIT is known to coincide with NEO).
To compare the simulated and measured Ar density, shown in figure 6(f ), recall that the CXRS diagnostic observes a single line, in this case Ar 16+ .This line is typically the most abundant at AUG temperatures, but lower ionization states are more populated in the edge, at lower T e , and the two higher ionization states dominate in the hotter core.This is why the measured Ar 16+ is hollow in the centre.This line has gradients and content inconsistent with the simulated total Ar density, but the correct comparison is against the simulated Ar 16+ density, which is provided in output by STRAHL.For Ar 16+ , the simulation matches the shape and content of the measured profile quite well, with a slight shift from the core to the edge which could be explained by the higher simulated core T e and differences in the Ar transport.
In figure 7 we expand on the simulation of radiated power densities for this discharge.2D tomographic reconstructions of the total radiated power density from bolometry were available for AUG #37041, and are shown in figure 7(a).An X-point radiator (XPR) [102,103], caused by the increased Ar seeding and located outside ρ tor > 0.85, can be identified.By comparison, not much radiation is seen in the core (although the uncertainty of the bolometry increases towards the centre due to lower line-of-sight coverage).The 2D SXR tomography is shown in figure 7(b), where we see that the experimental SXR emissions increase towards the centre but decrease again deep in the core, and there is a slight in-out poloidal asymmetry.
Figures 7(c) and (d) present results of the 1D ASTRA simulations in comparison to flux surface averages of (a) and (b), respectively.In figure 7(c), the Ar radiation is the dominant component and it peaks at ρ tor ≈ 0.8, where Ar is partially ionized but has a higher content than farther out radially, due to the edge gradient.W radiates only in the core.Note that a separate X-point contribution to the radiation is added.This is described in section 4.3.In total, the simulation obtains a high radiated power of 4.4 MW, in close agreement with the 4.5 ± 0.5 MW estimated by bolometry.Note that P LH and P sep are matched at 2.5 MW.The shape of the SXR radiation profile, in figure 7(d), is well reproduced from edge to ρ tor ≈ 0.25, where n w peaks and with it the modelled SXR emissions.

Model for radiation localized above the X-point
2D processes that produce a cold, dense region at the Xpoint (power balance between parallel heat conduction from the upstream and losses due to atomic processes [104]) cannot be modelled by ASTRA, which adopts a FSA approach to the transport equations.The formation of an XPR causes high radiation in a small edge region.Omitting this peripheral radiation leads to powers radiated inside the separatrix that are considerably lower than the measurements, while compensating with higher core radiation would trigger, by experience, a radiative collapse in the simulation.
We implemented a simple model to account for this peripheral radiation, estimating its magnitude as where L z is the total impurity cooling factor, n x and T x are the electron density and temperature of this dense, cold region above the X-point, c zx is the impurity concentration in this region and V x is its volume.Pressure balance along the field lines yields n x ≈ n u T u /T x , where the subscript 'u' denotes upstream values, here taken as the separatrix values in ASTRA.The key assumption of this model is therefore the value of T x , for which we take 7 eV as a typical value in the range of 1-10 eV obtained in AUG SOLPS-ITER simulations [105] and divertor Thomson scattering measurements [106].Two further assumptions are that the X-point impurity concentration is equal to its upstream value, i.e. the impurity density is compressed by the same amount as the electron density, and the value for the volume of the region (taken as a ring of circular cross section with radius R 0 ).The chosen parameters are illustrated in figure 8.We obtain 0.6 MW of additional radiation, added to the ASTRA simulation in a Gaussian profile with location and width in ρ tor of 0.97 and 0.10 respectively.

Impact of impurities on confinement
Impurity dilution has long been known to stabilize ITG turbulence [107][108][109].This is the explanation for the H-modelike confinement of this radiative L-mode provided in [94].
We investigate this further in this section by performing an additional simulation without impurities, such that n i = n e and Z eff = 1, but keeping the radiation profile fixed so the power balance is unchanged.TGLF is run with two kinetic species instead of five.In this case there is no dilution, the cause of ITG stabilization by impurities.Figure 9 shows the resulting simulations with and without impurities.The pressure profile p = n e T e + (n i + n B + n Ar + n w ) T i , whose volume integral corresponds to the stored thermal energy and ultimately characterizes confinement, is plotted in figure 9(a).The one with impurities is higher than the one without even though the main ion density of the latter is not diluted and the electron densities are equal, which means that the difference is caused by increased temperatures.Both edge pressures (ρ tor > 0.8) are practically identical, indicating that this is a core effect.The predicted ion and electron heat diffusivities, χ i and χ e , are shown in figure 9(b).Both are reduced from the pure to the impure plasma, particularly in 0.3 < ρ tor < 0.8.The result is a 33% increase in the predicted confinement, with H 98 = 1.06 and 0.80 with and without impurities.
To look deeper into the reduction of turbulent transport by impurities, we performed two standalone TGLF simulations using input data taken from radial averages in the yellow region in figure 9, where the pressure gradients differ most (R/L p of 25 and 17 with and without impurities).In particular, we build the inputs using the ASTRA simulation with impurities (red in figure 9) and then remove the impurities for the second standalone run, to have a more direct comparison of the role of impurities on the turbulence at the exact same background gradients.The corresponding turbulent spectra are shown in figure 10.The real eigenfrequency ω r , in figure 10(a), indicates turbulence propagation in the ion drift direction (negative ω r in the TGLF convention) at the larger length scales for both cases.At small scales, the simulation with impurities presents propagation in the electron drift direction but the one without impurities has parts of the spectrum on both signs of ω r .In figure 10(b), the growth rate γ shows three peaks, at large, intermediate and small scales.All three are reduced going from the pure to the impure plasma.Note that both ion and electron heat fluxes (Q i , Q e ) peak at larger length scales, around a normalized wavenumber of 0.3 characteristic of ITG turbulence, as shown in figure 10(c).Heat flux contributions from small scales are negligible.The inclusion of impurities reduces Q i and Q e by approximately a factor of 3.
The increased performance due to ITG stabilization by impurities begs the question of whether this effect would balance the deleterious fuel dilution in a fusion reactor.Fusion power grows with the square of the fuel density, P fus ∝ n 2 DT .In our simulation, the core fuel is diluted to 75% of the electron density.Given the 33% increase in confinement, this requires a scaling of P fus with at least the square of the confinement (which was found numerically for the confinement enhancement factor, P fus ∼ H 2. 15  98 , in [110]).However, the high confinement of this radiative L-mode is due to increased core temperatures and, as pointed out in [94], the D-T cross section decreases with T i at high enough temperatures.Still, this is an interesting scenario that is both ELM-free and has a high core radiated fraction (P rad /P aux ∼ 0.67), which are essential requirements for a safe power exhaust in a future reactor.

H-mode phase of this discharge
As previously mentioned, AUG #37041 had an initial H-mode phase before transitioning back to L-mode when the high Ar seeding was injected.Quasilinear transport models are typically not applied in full-radius H-mode simulations because they usually do not form strong pedestals.However, we saw in section 3 that suitable E r boundary conditions can result in pedestal-like structures in simulations with TGLF-SAT2.
As an exploratory endeavour, we performed a full-radius simulation of the H-mode phase of AUG #37041 (at 2.0-2.5 s), allowing E r to evolve freely up to the separatrix, i.e. 'floating' in equation (1).Note that using this choice in E r in the radiative L-mode phase would result in higher edge gradients than what was obtained in figure 6, overestimating the confinement in that case.In the experiment only the Ar puffing rate is different, leading to a reduced radiation and thus higher P sep > P LH .All other engineering parameters (P aux , I p , B φ , n e , are constant.For this simulation, we use the same B source as in the radiative L-mode phase, whereas both the and W sources reduced by a factor of 10, considering that in the experiment the Ar puff rate increased by around this much going from the H-mode to the radiative L-mode. Figure 11 presents the resulting simulated profiles.The main plasma profiles are of particular interest in this case, because we are applying TGLF-SAT2 up to the separatrix in an H-mode.We find surprisingly good agreement between the simulation and the experiment.Strong edge temperature and density gradients are obtained in figures 11(a) and (b), consistent with the experimental data.All three profiles of T e , T i and n e are within the error bars of the measurements across the full radius.In figure 11(c), Z eff also matches the IDZ estimate, which is about two thirds of the one in the radiative L-mode phase.
These H-mode-like profiles of the main plasma allow us to test the impact of the presence of a pedestal on the impurity density profile predictions.
There is a strong simulated edge n B gradient, as shown in figure 11(d), that unfortunately cannot be compared to measurements due to the lack of edge CXRS data in this phase (this is also evident in the T i measurements).The simulated profile is peaked, in contrast to the experimental profile which is flat at mid-radius and hollow in the core, and the B content is well reproduced.The W profile shown in figure 11(e), in contrast to the radiative L-mode phase, has a very steep edge gradient arising from the higher neoclassical pinch due to the n e pedestal and the lower edge turbulent diffusivity.It is then hollow in the region 0.6 < ρ tor < 0.8 (where n e is flat while T i is peaked, so there is temperature screening), and it finally peaks in the core, although less than in the radiative L-mode phase.The magnitude of the effect of adding collisions with Ar are smaller than in figure 6(e) because the Ar content is lower in this phase, with a qualitatively similar increase in the core peaking.Both measured and simulated Ar contents, shown in figure 11(f ), are about 1/6 of the ones in the radiative L-mode phase, even though the seeding rate was reduced by a factor of 10, because Ar also presents a stronger edge gradient that increases the core density from its separatrix value.
In figure 12(a) similar analysis of the radiation as in figure 7 is presented.Figures 12(a) and (b) show the tomographic inversions of the measured total radiation density and its SXR share of the spectrum.Note the lack of the XPR in this phase of the discharge.The contribution of W to the radiation is twice as in the L-mode phase even though the source is reduced by a factor of 10, as shown in figure 12(c), because the W density is more distributed in the edge.Ar radiation is reduced by a factor of 5.The total 2.3±0.5 MW estimated from bolometry is well matched by the simulated 2.2 MW.In figure 12(d) we see that Ar still contributes significantly to the SXR emissions.

Discussion
Full-radius integrated modelling of ASDEX Upgrade L-mode experiments has been performed, with particular emphasis on the transport of multiple impurity species and their effects on the main plasma through radiative cooling and fuel dilution.This has been achieved thanks to the development of a modelling framework that integrates TGLF-SAT2, FACIT and STRAHL in ASTRA, applying the quasilinear turbulent transport model from axis to separatrix to obtain impurity transport coefficients for the first time.
A set of six L-modes, including intrinsic boron and tungsten impurities in the plasma, was simulated to test the predictions of the workflow against experimental profiles under differing plasma currents and heating power mixtures.Important effects, such as the increase of confinement with plasma current and the control of core tungsten accumulation with central ECRH, were reproduced in the modelling.
The choice of edge sources, which are currently ad-hoc inputs used to approach expected radiated powers and effective charges, is an evident aspect to improve upon.However, an accurate description of sources requires more sophisticated SOL modelling capabilities than STRAHL's 1D model.While initial efforts in integrated core-pedestal-SOL modelling have recently begun, the computational cost of this type of simulations remains high.Likewise, a model to accurately predict the neutral density profile at the edge is necessary in order to include missing physics like the effect of CX reactions on the peripheral impurity distribution and radiation.The ionization electron source from partially stripped impurities could also be included by updating the STRAHL coupling to ASTRA.A more sophisticated treatment of MHD transients is another element to improve upon.Moreover, momentum transport modelling is a missing element that would provide the toroidal rotation profile self-consistently.This is a high priority development for future work, given the strong effects of rotation on the neoclassical and turbulent transport of heavy impurities.
A radiative L-mode scenario with high confinement and no ELM activity has been modelled in close agreement to multiple diagnostics measuring the main plasma, impurity densities and radiation.The effect of ITG stabilization due to impurity dilution of the fuel has been confirmed to play a significant role on the enhanced confinement of this experiment.A simple XPR model has been applied to account for highly localized radiation at the edge.While introduced here retroactively (that is, knowing the XPR is present by looking at the bolometry), the XPR access condition of [104] could be implemented in subsequent integrated modelling applications.Furthermore, the E×B shearing treatment of [39] proved sufficient to reproduce experimental profiles in the H-mode phase of this particular discharge, giving hope to further applications.
The radiative L-mode scenario could be further developed, in a predict-first approach with simultaneous integrated modelling efforts like the one performed here.More efficient radiators [4], active XPR control [111] and the new AUG upper divertor [112] (starting late 2024) could be leveraged in future experiments.
There is a need to further validate the treatment of multiimpurity collisions in FACIT against NEO, possibly extending parts of the model to better reproduce the drift-kinetic results.A precise description of multi-impurity effects with toroidal rotation still requires the use of NEO, which carries significant consequences on the computational time of the integrated modelling framework.
A clear objective for future work is to extend the applicability of the framework to full-radius simulations of seeded H-mode plasmas.This could be done by including a description of impurities in IMEP [12,13].Three key elements to be addressed, moving from L-mode to H-mode simulations, are the estimation of turbulent impurity transport coefficients in the pedestal, the impact of impurities on the SOL model that provides boundary conditions at the separatrix in IMEP, and effects of high impurity content and radiation on the pedestal stability.

Figure 1 .
Figure 1.Schematic representation of the modelling framework including impurities.The blocks outside of the dashed lines represent external user inputs, whereas everything inside them is self-consistently integrated in ASTRA.

Figure 2 .
Figure 2. ASTRA simulation of AUG #39255 at 2.5-3.0 s (high current, mixed heating).(a) Electron (red) and main ion (blue) temperatures, with measurements in discrete points and simulations in solid lines.(b) Electron and main ion densities, colour-coded as in (a).(c) Effective charge and safety factor simulated profiles in the solid pink and brown lines.(d) Total simulated boron (solid) and B 5+ line (dashed) densities in green.(e) Tungsten density simulation in the solid cyan line.Measurements in (c)-(e) are shown in black dotted lines with uncertainty bands in grey.The scaling factors of the experimental data in (d) and (e) are discussed in the text.(f ) Radiated power density calculated by STRAHL (red) with contributions of W (cyan) and B (green); the core radiated power, its components and measured value are in the label.

Figure 4 .
Figure 4. (a) Stored thermal energy of the plasma as a function of the plasma current.(b) Boron density gradient as a function of the electron density gradient, both at mid radius (averaged over 0.45 < ρtor < 0.55).(c) Tungsten peaking (on-to off-axis nw ratio) as a function of the ECRH power fraction.Solid symbols are the simulation results, while open symbols are the experimental data.The different colours represent different heating mixes, and the different markers correspond to different plasma currents.

Figure 5 .
Figure 5. Choice of the boundary condition of the simulated radial electric field in AUG #39323 at 5.0-6.0 s.(a) Edge Er profile at different values of ρ min .The inset figure is a zoom of the y-axis.(b) Resulting edge main ion pressure profile, with the experimental data in black and its uncertainty band in grey.

Figure 7 .
Figure 7. Radiated power density measurements and simulations for AUG #37041 at 5.0-5.5 s.(a) 2D tomography of the bolometry.The solid red contour is the separatrix and the dashed red contour marks ρtor ≈ 0.85.(b) SXR measurements.Here the dashed red line is the surface at Te ≈ 1 keV.(c) Radiated power profiles: FSA of bolometry from (a) in black, simulated radiation in red with its components by W (cyan), Ar (purple), B (green) and the X-point radiation (yellow, see section 4.3).The volume integral of each contribution is given in the label.Vertical lines correspond to the contours in (a).(d) SXR power profiles: FSA of measurements from (b) in black, with reconstruction from the modelling using temperatures and densities from figure 6 and SXR cooling factors.

Figure 8 .
Figure 8.(a) Bolometry data for X-point radiation, with cross sections of possible volumes to be used in equation (9) shown in circles.(b) Modelled X-point radiated power as a function of the X-point temperature for the three volumes shown in (a); the selected parameters are highlighted in the red circle.

Figure 9 .
Figure 9. (a) Simulated pressure profile for AUG #37041 at 5.0-5.5 s without (blue) and with impurities (red).(b) Main ion (solid) and electron (dotted) heat diffusivities.The yellow region is where the turbulent spectra will be analysed in figure 10.

Figure 10 .
Figure 10.TGLF spectra for AUG #37041 at 5.0-5.5 s at 0.70 < ρtor < 0.78.(a) Real eigenfrequency of the dominant mode.Negative and positive values represent the ion and electron drift directions.(b) Growth rate of dominant mode.Both variables in (a) and (b) are normalized to the sound speed over the minor radius and divided by the wavenumber.(c) Electron and main ion heat fluxes in gyro-Bohm units.