SOLPS-ITER modeling with activated drifts for a snowflake divertor in ASDEX Upgrade

We report on the first SOLPS-ITER simulations of a low-field side snowflake minus (LFS SF−) divertor configuration with drifts fully activated in ASDEX Upgrade. Compared to a reference case without drifts, the simulation in normal toroidal magnetic field configuration (B × ∇B points to the primary X-point) shows a larger low-field-side/high-field-side asymmetry, an enhanced radial cross field transport, as well as a flux redistribution between the primary and secondary strike points. Although small compared to the total input power, power is found even on a strike point magnetically disconnected from the outer mid-plane, which is hard to explain by purely diffusive transport.


Introduction
In a magnetically confined fusion device, the material limits of the plasma facing component can impose serious restrictions on the plasma operation window. ITER will be operated in a partial detachment regime to avoid intolerable power and particle fluxes to the divertor target. For future fusion reactors with even higher power, alternative divertor configurations may be necessary to further reduce the target load. For this reason, various alternative configurations, such as the snowflake divertor [1], the X-divertor [2], and the super X-divertor [3], were designed and experimentally tested in recent years. In order to study the physics in these configurations in detail and predict their performance in future devices, numerical simulations were carried out applying 2D and 3D codes, such as SOLPS [4], SOLPS-ITER [5], EDGE2D [6], UEDGE [7], SOLEDGE2D [8] and EMC3-EIRENE [9].
For the ASDEX Upgrade tokamak (AUG), a modification of the upper divertor is currently prepared in order to study various alternative divertor configurations experimentally in a machine with a high heating power compared to its size [10]. Recently, the SOLPS5.0 code package was successfully applied to simulate a low-field side snowflake minus (LFS SF − ) configuration based on the new upper divertor geometry [11]. It was shown that divertor detachment was reached either with lower upstream density or with lower impurity seeding rates in the snowflake configuration compared to a single null (SN) reference with the same input parameters. However, drifts were not included in this work, which were expected to contribute to a detachment asymmetry between the inner and outer divertor and a modification of the cross field transport as found in experiments and simulations [12,13]. So in order to make realistic predictions, it is important to activate drifts and study their effects in such simulations. One of the main questions related to the benefits of a snowflake configuration are concerning its capability to share heat and particle fluxes among its strike points. Drift effects, e.g. increased radial transport and modified parallel impurity transport, as well as predicted enhanced ion losses, are expected to play an important role [14]. In EMC3-EIRENE simulations for the snowflake configuration in TCV, the modeling without drifts [15,16] underestimated the heat flux at the secondary strike point by about one order of magnitude compared to experimental measurements [17]. In a subsequent work [18], a qualitative explanation for the activation of the secondary strike point was found by computing drift terms on the plasma background given by EMC3-EIRENE, however, a quantitative and self-consistent simulation is still missing.
It is well known that SOLPS-ITER simulations with drifts can be numerically very expensive. Fortunately, a speed-up method [19] developed recently made it possible to achieve convergence of the simulation within a few months. In this paper, we present the first simulation results for a snowflake configuration with drifts fully activated in SOLPS-ITER.

Simulation setup
While SOLPS5.0 was applied in a previous study of a snowflake configuration without drifts [11], here we apply the SOLPS-ITER 4 code package which contains a more complete treatment of the drifts and currents. It consists of the multifluid transport code B2.5 [20] and the Monte Carlo code EIRENE [21] simulating the neutral particle transport. The simulation mesh for the plasma is very similar to the one used in [11], while a triangle mesh for kinetic neutrals required by SOLPS-ITER is added. The plasma is assumed to be pure deuterium and the anomalous radial diffusion coefficients (D n AN , c e AN and c i AN ) to be spatially constant. A simulation with impurities and a radial transport barrier requires more computational time to converge and is foreseen in the future.
In this work, E×B and diamagnetic drifts, as well as currents caused by viscosity, ion-neutral friction and inertia are all fully activated, and the equations are chosen according to [22] where a reasonable trade-off was made between physical accuracy and numerical stability.
The equations are solved on a 2D grid, where x and y represent the poloidal and radial grid coordinates, respectively. The ion particle continuity equation is given by ⎛ where n and S are the ion density and source rate, h x , h y and h z =2πR are cell lengths in poloidal, radial and toroidal directions, respectively, and = g h h h x y z . The particle flux in the poloidal (x) and radial (y) directions can be written as˜( The first terms on the rhs are given by where D n is the particle diffusivity and b x =B x /B is the ratio of the poloidal component to the total magnetic field. Considering the classical diffusivity to be negligible compared to the anomalous one (D n AN ), and applying a numerical correction with the convective velocity disscussed in [22], D n is set as In the poloidal direction, besides the diffusive flux, Γ 0,x also contains a convective flux from the poloidal projection of the parallel velocity and the Rhie and Chow upwind correction velocity [23] ⎛ where m is the ion mass, and p=n(T e +T i ) the thermal pressure.
The second terms are the fluxes driven by the E×B drift where f is the plasma potential. The particle fluxes related to diamagnetic drifts are given by . 8 Since equation (1) only depends on the divergence of these terms and not on the actual value itself, in the numerical procedure the diamagnetic drift is replaced by an effective drift (equation (9) in [22]) that has approximately the same divergence. This is done in order to improve the numerical stability, while it preserves the result. When we discuss the role of the diamagnetic drift later on, however, we compute the fluxes according to equation (8) by post-processing the converged solution.
Currents caused by anomalous transport ( j AN ), inertia ( j in ) and viscosity ( || j vis ,j vis⊥ and j visq ) are assumed to contribute to the particle flux in radial direction only While the poloidal component is set to zero, Γ current,x =0, for reasons of numerical stability.

Simulation results
The LFS SF − configuration simulated in this work is shown in figure 1. A secondary X-point is located in the low-field side scrape-off layer (SOL) dividing it into near and far-SOL regions connected to the areas around OT1 and OT2, respectively. In addition to that the two remote regions R1 and R2 are connected to a third strike point at OT3. The x (poloidal) and y (radial) directions in the different regions of the grid are marked by the black and white arrows, respectively. The arrows point into the positive direction, along which the cell index increases. The deuterium ion density and the electron and ion temperatures are fixed at the innermost simulation boundary at 2×10 19 m −3 and 200 eV, respectively. The boundary condition for D + at the outermost radial boundary is set by a leakage flux Γ leakage =0.001nc s , where n and c s are the local density and sound speed, respectively. The recycling coefficient at the cryo pump surface (see figure 1) was set to 0.9, corresponding to a pumping speed about 50 m 3 s −1 .
The deuterium atomic and molecular processes including ionization, dissociation, charge exchange, elastic collisions and recombination were simulated by EIRENE. The number of test particles used in the Monte Carlo procedure was 57 000, and an averaging scheme [24] was applied to decrease statistical errors.  figure 2 with a temperature decay length about 10 mm in the near SOL. This decay length is a common value in AUG attached discharges [25]. The toroidal magnetic field points into the paper leading to an upward B×∇B drift towards the upper divertor. The convergence criteria used in this work are based on those described in [5], including the steadiness of density, temperature, parallel velocity, potential and total recombination rate on a typical radial diffusion time scale (0.2 s). In addition to that, the error of the global particle flux balance is required to be less than 1% and that of the energy balance of order 1% (see the last rows in tables 1 and 2). The largest residuals of the power transport are located in regions where the grid resolution is coarse (e.g. the inner divertor and the region around the primary X-point).
In order to show the effect of drifts clearly, here we compare the case with drifts activated to a case without drifts but with the same input parameters otherwise. The density, temperature, pressure and potential profiles at the outer midplane in these two cases are shown in figure 2. Compared to the non-drift case with the same transport coefficient = D 0.5 n AN m 2 s −1 (green dashed line in figure 2(a)), the density profile in the drift case (blue line in figure 2(a)) is flatter due to the drift effect on the radial transport. The different upstream density then makes the comparison at the target difficult. For this reason, the particle transport coefficient in the non-drift case is increased to = D n AN 1.0 m 2 s −1 to achieve similar upstream density, temperature and pressure profiles (black dashed line in figure 2), while c e AN and c i AN are kept at 1.0 m 2 s −1 . In the case with drifts, the potential radially increases inside the separatrix and decreases in the far-SOL, which leads to a sheared radial electric field (pink line in figure 2(d)). The power and particle balances in the simulations with and without drifts are shown in tables 1 and 2. It is found that the power heating the innermost radial boundary (Q CORE ) and leaving the outermost boundary (Q WALL ) are both slightly higher in the non-drifts case, while the power to the targets (Q TAR ) are similar in these two cases. The difference in the particle flux into the innermost boundary is negligible compared to the ionization source. Figure 3 shows the profile of the poloidal ion particle flux caused by drifts at the outer mid-plane. The poloidal E×B and diamagnetic fluxes are shown by blue and green dashed lines, and the total drift flux (Γ drift,x =Γ E×B,x +Γ dia,x ) is plotted by the black solid line. Near the primary separatrix (r u =0), Γ E×B,x and Γ dia,x have opposite signs and tend to cancel each other. In the far-SOL, Γ E×B,x and Γ dia,x are pointing into the same direction and lead to an anticlockwise  poloidal drift flux, i.e. poloidally from the high-field side to the low-field side along the poloidal coordinate in the SOL.
To investigate how drifts impact the divertor target behavior, we show the distribution of the total drift particle flux in the divertor region in figure 4, and the target profiles of the parallel energy flux, density, temperature, total pressure and current in figure 5. The target profiles are compared with those in the case without drifts (shown by black dashed lines). Figures 4(a) and (b) show the poloidal and radial projections of the total drift flux. In the SOL, the poloidal drift flux is mainly transporting particles from the high-field side to the low-field side where it is split between OT1 and OT2. Since the particle fluxes are taken into account as a convective term in the energy flux expression [22], it is expected that the drifts have an impact on the energy flux. This is indeed seen in figures 5(a) and (f) which show a lower energy flux at the inner target and higher ones at OT1 and OT2 near the strike points in the drift case. The parallel energy fluxes at the targets were estimated by  » q q b , where E pot is the average potential energy released by a deuterium ion recombining to a molecule at the target surface. The expression includes non-parallel contributions from drifts. The poloidal electron heat flux q e,x at the target corresponds to the sheath current-voltage characteristics and the full expression can be found in equation (35) in [22] and the appendix C.7.4 in [26]. Near the inner target, due to the large poloidal temperature gradient, a strong radial outward drift flux is found (see figure 4(b)). It pushes the high density region near the inner strike point towards the far-SOL and cools down the outer part of the inner target, which can be seen in figures 5(b) and (c). A region of high density and low temperature in the high-field side far-SOL is seen in the 2D plot in figures 6(a) and (b), in contrast to the non-drift case (figures 6(c) and (d)).
In the high-field side SOL of the snowflake configuration, the E×B and diamagnetic drifts act quite similar to those in a SN geometry, while in the low-field side SOL, the situation becomes more complicated. The drifts not only simply enhance the total energy flux to the outer targets, but also have an impact on the distribution of the energy to the three outer targets. As shown in figure 5(f), due to the radial inward drift in the low-field side SOL (see figure 4(b)), more energy is transported to OT1 while the outer part of OT2 receives slightly less energy compared to the non-drift case. The region near OT3 is zoomed in figure 4(c), and the arrows point into the direction of the total drift flux. Due to the radial flux across the OT2 divertor leg and the poloidal flux from OT2 to OT3 in the region R2, the peak value of the energy flux at OT3 increases by a factor of 4 and becomes comparable to that at OT2. Such an activation of OT3 was observed experimentally in TCV [15,17] and qualitatively explained by computing drift terms on the plasma background given by EMC3-EIRENE simulations [18]. In the region R1, the poloidal drift flux from OT3 to OT1 partly counteracts the increase of energy flux to OT3 and contributes to a small peak in the outer part of the OT1 energy flux profile (s ot ≈22 mm). In the case with drifts, the residual of the power transport summed over the cells in the regions R1 and R2 is 0.08 kW, which is much lower than the integral power to OT3 (9 kW). In the SOL near the targets and just outside the separatrices, reversed drift fluxes are found (see figures 4(a) and (c)), which could contribute to multiple peaks in the energy flux profiles. Multi-peak target profiles were also observed in SN configurations in JET [27] and in snowflake configurations in TCV [17]. EDGE2D [28] and SOLPS5.2 [29] simulations showed that the drifts played an important role to explain this phenomenon.
The change of the temperature at the outer targets is the same as the change of energy fluxes, i.e. the temperature at OT1 and OT3 increase and the one at the outer part of OT2 slightly decreases (see figure 5(h)). It should be noticed that although the maximum temperature at OT1 (about 30 eV) is increased by drifts, it is still lower than the one (about 45 eV) in a SN reference simulation without drifts as shown by the cyan dotted lines in figure 5(h). For the sake of brevity the SN reference simulation is not shown here explicitly. However, a plot of the computational grid can be found in figure 1 of [11]. The transport coefficients and boundary conditions used for the SN reference simulation are the same as those in the snowflake non-drift case. This results in nearly identical upstream profiles in these two cases (see cyan dotted curves in figure 2). The power and particle balances for the SN reference are given in tables 1 and 2 in the last column. Similar to the snowflake case, one can expect that switching on the drifts Table 1. Power balance in the simulations with and without drifts. Q CORE is the power heating the innermost radial boundary and Q WALL is the one leaving the plasma through the outermost radial boundary of the regions FARSOL, PFR, R1 and R2 (see also figure 1). Q TAR is the total power to the divertor targets IT, OT1, OT2 and OT3. Q rad is the total radiation power. The balance for a single null (SN) reference simulation without drifts is shown in the last column.

Power
SF with drifts SF without drifts SN without drifts would also increase the maximum temperature at the outer target in the SN case, as discussed in [30]. Without considering the recycling, particles transported poloidally to the outer divertor by drifts could increase the density there. However, in the situation in this paper, the density near the divertor target is dominated by the ionization of recycled particles. In the drift case, the higher power and temperature near OT1 lead to a lower collisionality (n~n T 3 2 ), the strong reduction of the ionization source results in a much lower density (see figure 5(g)). The radial drift from the low-field side SOL across the separatrix near the primary X-point and the OT1 divertor leg to the high-field side could also contribute slightly to the reduction of the density, as well as the pressure (see figure 5(i)) at OT1. Figure 5(e) and (j) show the total current perpendicular to the inner and outer targets, respectively. In the SOL, the current is dominated by the poloidal projection of the parallel thermal current and travels from the hotter divertor (OT1 and OT2) to the colder one (IT). In the case with drifts, reversed currents are found in the narrow region near the strike points.

Conclusion and outlook
Converged SOLPS-ITER simulations with drifts fully activated are carried out in a LFS SF − divertor configuration for the first time. The configuration corresponds to that of the future upper divertor of ASDEX Upgrade in normal toroidal field direction (B×∇B points to the primary X-point).
Comparing with a reference case without drifts but with similar upstream profiles, the simulation with drifts shows a larger low-field-side/high-field-side asymmetry and an enhanced cross field transport especially in the inner divertor. In addition, a redistribution of particles and energy among the primary and secondary targets is found in the outer divertor. With drifts, the target OT3 magnetically disconnected from the outer mid-plane is activated reaching a maximum parallel energy flux density by a factor of 4 higher than that in the non-drift case, although the integral power to OT3 is still small compared to the total input power. The activation of a secondary strike point like OT3 was one of the motivations to study snowflake configurations [1] and was  observed in experiments [15,17] and non self-consistent simulations [18] in TCV. Given that the integral energy flux to OT3 remains small, the relevance for a reactor still needs to be tested by future simulations studying how this effect scales with β and/or machine size.
In this work, the plasma was assumed to be pure deuterium, while the impurity radiation and transport with drifts can also be important for evaluating the power exhaust capacity of such a snowflake configuration. This is foreseen as the future work.