PSFC/JA-20-59 Simulations of divertor heat flux width using transport code with cross-field drifts under the BOUT++ framework

The authors would like to thank T. D. Rognlien, B. Chen, J. G. Chen, Zeyu Li, T. Y. Xia, Y. M. Wang, and T. F. Tang for useful discussions and A. Q. Kuang and D. Brunner for providing heat flux data for discharge (Grant No. 1160729008). This work was supported by the National Key R&D Program of China (Grant Nos. 2017YFE0301206 and 2017YFE0300402) and the National Natural Science Foundation of China under Grant Nos. 11675037 and 11575039. This work was also performed under the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract No. DE-AC52-07NA27344. Further support comes from US DoE (Award No. DE-SC0014264). The LLNL work was supported by the U.S. DoE, Office of Science, Office of Fusion Energy Sciences. Reproduction, translation, publication, use and disposal, in whole or in part, by or for the United States government is permitted.


Introduction
The high confinement mode (H-mode) features a steep pressure gradient at the plasma edge, called the pedestal region. The H-mode is the baseline operation regime for the high-performance requirements of future steady state tokamaks such as ITER 1, 2 . Diffusive and convective transport of plasma from the strong gradient region of the edge pedestal into the scrape-off-layer (SOL) can impact power exhaust, which will be a critical issue for future magnetic confinement fusion devices. From power balance in steady state, the power across the separatrix must equal to the input power, plus fusion alpha particle power, minus core radiated power. As the alpha power in burning plasmas is large, this can lead to large power fluxes on the divertor or other material surfaces, leading to damage of the plasma-facing components (PFCs) 2 . The width of divertor heat flux footprint, or so called scrape-off layer (SOL) width ! , and the peak of the heat-flux density on divertor targets are critical parameters for characterizing the particle and power exhaust 3 . The size and scaling of the heat-flux width have been the subject of recent experimental investigations and an inverse scaling dependence of the heat flux width on plasma current Ip is observed from a data set combining multiple machines [4][5][6][7][8] , including NSTX, ASDEX Upgrade, Alcator C-Mod, JET, and DIII-D H-mode discharges. A database from these experiments gave rise to a scaling law by Eich, in which ! is nearly inversely proportional to the poloidal magnetic field at the outer midplane, with weak dependencies on any other key SOL parameters 9 . A similar dependence is also found in the EAST tokamak 10 . From this scaling, a very narrow SOL inter-ELM heat flux width, !~1 , is projected for the ITER baseline burning plasma scenario (Ip=15MA). The divertor target power loads with such narrow width would in reality be impossible to handle, which raises a serious concern 2, 11 . Even though the divertor detachment is the most promising means for mitigating target heat fluxes, the narrow lambda_q|| will make the plasma harder to detach.[ Kukushkin A.S., Pacher H.D., Pacher G.W., Kotov V., Pitts R.A. and Reiter D. 2013 J. Nucl. Mater. 438]. To date, no consensus has been reached for the magnitude of the heat flux width for ITER and future reactors. Additional recent results on C-Mod have shown consistency with the Eich scaling up to ITER-relevant Bp of 1.3T 12 . However, gyro-kinetic simulations and turbulence simulations at the ITER scale indicate !~5 − 6 [13][14][15] . No matter which prediction is right, there is no doubt that the heat flux width definitely needs to be broadened by any applicable techniques.
For a given input power, the peak of heat flux is inversely proportional to a critically important parameter -the scrape-off-layer heat flux width. It is generally thought that ! is determined by the competition between perpendicular transport governed by drifts and turbulence and parallel transport in the SOL. In order to explain the physics behind the inverse scaling of ! with Ip, Goldston proposed a heuristic drift-based (HD) model, which is consistent with the international multi-machine scaling. The hypothesis in the HD model is that magnetic drifts balance the nearly sonic parallel flow in setting the ion particle flux width and cross-field transport at the separatrix is dominated by classical drifts 16 . In the HD model, the power across the separatrix is empirical and ExB drift is not considered. Recently, electrostatic kinetic simulations with the PIC code XGC reproduced the trends in the Eich scaling and HD model for existing tokamaks 13,17,18 . In addition, the turbulence code (6-field two-fluid module) [19][20][21][22] under the BOUT++ framework has been developed to study the edge physics in tokamaks, including ELMs and the divertor heat fluxes during the ELMy Hmode. The role of turbulence in setting the electron heat flux width has been studied for the inter-ELM phase of DIII-D 23 , EAST ELMing H-mode discharges 24 and C-Mod EDA H-mode discharges 25,26 . The turbulence simulations of C-Mod discharges further demonstrate that reducing the radial electric field Er, which is calculated from the transport module under the BOUT++ framework without the cross-field drifts, can effectively enhance the radial turbulence transport, thus increasing the heat flux width 26 .
The BOUT++ turbulence code mainly focuses on fast turbulence dynamics and scale of up to a few thousand Alfvén times. Due to slow ion parallel dynamics, even though the turbulence simulations are able to reach a steady state solution at the midplane for both ions and electrons, and electron divertor heat flux is able to reach a well saturated state, the ion heat flux on the divertor targets has not been easily characterized for the BOUT++ turbulence simulations. For quick scoping studies, a new transport module, including all drifts, which mainly focuses on 2D plasma evolution to steady state on a longer transport time scale, is developed under the BOUT++ framework 27,28 to get a steady state solution of divertor heat flux for both electrons and ions. In this work, we will use the new transport code and the same 4 discharges from C-Mod as in Refs. 25 and 26 (1) to study the relative role of cross-field ExB and magnetic drifts vs. turbulent transport in setting the heat flux width and (2) to demonstrate the transition from drift to turbulence dominant regime. In addition, we will make the comparison of electric field between simulations and experimental measurements with or without magnetic drifts. These studies feature the following significant new results.
1) Compared simulated electric field Er with experimental measurements, as shown in Fig.5. 2) Compared simulated divertor heat flux width with measurements as shown in while the ExB drift decreases the heat flux width by 10%~25%, which gives improved agreement with experiment relative to Goldston's model, as shown in Fig.11. 6) Identified two distinct regimes via a turbulence diffusivity scan ( " ): a drift dominant regime when " is small and a turbulence dominant regime when " is large. Magnetic drifts give a lower limit of the width ! , similar to Goldston, as shown in Fig 12. These studies suggest that both drifts and turbulent transport compete to determine the divertor heat flux width in C-Mod EDA H-mode discharges.
The rest of the paper is organized as follows. Section 2 introduces the transport physics model with all drifts and the setup of the heat flux width simulations. The steady state Er has been compared with experimental measurements of a C-Mod discharge in section 3. The simulation results of heat flux for C-Mod EDA H-mode discharges are described in section 4. Importantly, the heat flux width calculated from the transport code (trans-Er module) has been also compared with the turbulence code (6-field module) results under the BOUT++ framework, experimental measurements and the Goldston HD model. The competition between drifts and turbulence in setting the divertor heat flux is shown in section 5. The last section summarizes the results.

Physical model and setup of simulations
In this paper, the transport model (trans-Er) with all drifts under the BOUT++ framework has been used to simulate the plasma transport and power deposition in C-Mod discharges 29 . Starting from the Braginskii equations 19,30,31 , the evolving variables are the vorticity ( ), electrostatic potential ( ), ion density ( # ), parallel ion velocity ( ∥# ), and electron and ion temperature ( ! and " respectively). The set of transport equations have been given in Ref 29 . For completeness, we write the equations here again: Here the ExB velocity is *0, = 5 , ! 2 × ∇ , and the ion diamagnetic velocity is Where the kinetic effects are introduced through C and C = 1.0 is used in the simulation. Here is the safety factor and 2 is the major radius. In the ion density and ion temperature equations, ExB drift and ion diamagnetic drift are considered. The ExB and electron diamagnetic drifts are included in the electron equation. The ion polarization velocity is a small correction to the ExB and diamagnetic velocities for transport problems, which is neglected in Eqs. (1)(2)(3)(4), but is kept in Eq.(5) of vorticity equation. This transport model with all drifts and radial electric field calculation are described explicitly as in reference 29 . It is worth noting that in the following simulations, the effective radial velocity from particle diffusion is included in the velocity, temperature and vorticity equations as part of convection, while the source terms associated with neutral particles are not included, such as ionization, recombination, charge exchange and radiation. The plasma fluid equations are discretized in a 2D orthogonal flux-surface mesh using a finite-difference algorithm.
Using a method similar to that in reference 29 , for a steady state solution without volume sources, drifts and flow, the ion continuity equation (1) becomes, (10) Given the core boundary particle flux and experimentally measured density profile inside the separatrix, the particle diffusivity D " can be obtained by D " Here the core boundary flux is determined by nominal boundary particle diffusivity " = 1 4 / and the experimentally measured density profile at the core boundary.
With the effective radial velocity from particle diffusion in the temperature equations, for a steady state solution without volume sources, drifts and flow, the electron and ion temperature equations (Eq. (3) and (4)) become, The heat transport coefficients "# and "( can be calculated by: E0 | 0G2 is given. In our simulations, at the core boundary x=0, the gradients are calculated from experimental profiles and a nominal thermal diffusivity "#,( =1m 2 /s is given. The sheath boundary conditions (SBC) are applied on the divertor targets surface which we assume that the ion parallel velocity matches the sound speed: , where # represents ion mass. The same SBC are used as reference 29 , which the total parallel current through the sheath is described as: The sheath potential can be derived as: In our simulations, the cross-field drifts (ExB and magnetic drifts) are included in the plasma transport equations, which are important in calculating the radial electric field (Er) and the parallel heat flux on the divertor targets 29 . To compare the Er calculation from our transport code with experimental results, a C-Mod EDA H-mode discharge (#1080321020) with lower single null divertor geometry is simulated. Fig. 1 shows the simulation domain which includes the plasma edge, the SOL region and private flux region. In the radial direction, the normalized poloidal flux is from ψ = 0.9 to ψ = 1.05, where ψ = 1.0 is the location of magnetic separatrix. The grid generated from the equilibrium file (EFIT g-file) and the spatial resolution of the mesh is 260 grid points in the radial direction and 64 grid points in the poloidal direction. Section 3 shows the results of this comparison.  Table 1: Key plasma parameters of the C-Mod discharges for 4 different plasma currents. The first column is the shot number. The second column is the time slice that the plasma profiles are used for the simulations. The third column is the toroidal magnetic field. The forth column is the poloidal separatrix magnetic field at outer midplane. The fifth column is the separatrix temperature as determined using unshifted Thomson scattering profiles and the standard EFIT magnetics reconstruction. The sixth column is the plasma current.
Next, in order to explore the physics behind the Eich scaling of the heat flux width ! versus Ip, such as drifts versus turbulent transport, a set of four C-Mod EDA Hmode discharges with lower single null divertor geometry are simulated. Table 1 shows the detailed key plasma parameters for 4 different currents with shot numbers #1100223012, #1100212023, #1100303017 and #1160729008. The simulation domain is the same as for the discharge in Fig. 1. The initial plasma profiles inside the separatrix of the four shots used in these simulations are fitted onto a radial coordinate of normalized poloidal magnetic flux from the experimental data with a modified tanh function. As an example, the squares in Fig. 2 are the experiment data points and the red curves are the fitting profiles of electron density and temperature for Ip=1.0MA case. in the simulations, outside the separatrix ψ > 1.0, the initial plasma profiles are assumed to decrease linearly from separatrix to outer boundary. The transport coefficients are calculated based on the initial experimental profiles. The flux-limited parallel thermal transport is implemented in the simulations. The initial profiles for four different currents in the simulations are shown in Fig. 3. The pedestal pressure profiles show a clear trend for these four shots as shown in Fig. 3(a). With the plasma current Ip increasing, the pedestal height increases or the scale length decreases, as seen in previous empirical results 32-34 .  One should note that the sheath boundary conditions (SBCs) are implemented at the divertor targets in the simulations as in reference 29 . The effects of cross-field drifts on the heat flux width are investigated by turning on or off the ExB drift and magnetic drift terms. In general, from the trans-Er module, we can derive: (1) the radial electric field; (2) radial and poloidal plasma profiles in the SOL; (3) radial and parallel heat fluxes profiles; (4) divertor heat flux magnitude and widths. The quantities that can be compared with experiments are the Er profile, the divertor heat flux magnitude and widths.

Er comparison with experiment
The onset of H-mode is characterized by a sharp decrease in edge fluctuations and the formation of edge transport barriers (pedestals) in both temperature and density profiles. These barriers lead to significant increases in core confinement as a result of temperature and density profile stiffness. It is commonly accepted that ExB velocity shear is responsible for the suppression of edge turbulence, which reduces the losses of both energy and particles across magnetic field. Classical particle motion from ExB and magnetic drifts is believed to be important for understanding tokamak edge/scrape-offlayer (SOL) transport. In this section, we perform transport calculations using the trans-Er module 29 under the BOUT++ framework. We also do the comparison of radial electric field with the experimental measurement for a C-Mod EDA H-mode discharge (#1080321020) with lower single null divertor configuration. The major parameters for this EDA H-mode discharge are: the toroidal magnetic field . = 5.4 T, the separatrix temperature !,0!1 = 63.5 eV, the separatrix poloidal magnetic field at outer midplane 1,234 = 0.66 T and the plasma current 1 = 0.81 MA. Figure 4 shows the initial density and temperature profiles we used in the simulations. For all C-Mod cases in this paper, ion temperature is assumed to be the same as electron temperature (Ti=Te) due to high collisionality. For the specific discharge considered in this Section (#1080321020), measurements from pedestal charge exchange recombination spectroscopy (CXRS) confirm the equivalence of Te and Ti, within experimental errors. The profiles in Fig. 4 are fits to experimentally measured data inside the separatrix, while in the SOL, they follow a constant absolute value decay from separatrix. The radial electric field is calculated by solving the potential (Eq. (7)) and vorticity (Eq.(5)) equations for two cases with and without cross-field drifts. In the simulation, #," is the anomalous perpendicular viscosity and #,∥ is the classical parallel viscosity. The effect of #," and #,∥ on the linear growth rate of peeling-ballooning model 35 and the steady state solution of radial electric field 36 has been discussed in the ELM simulation and Er calculation under the BOUT++ framework. The constant value of #," = 5 4 / and #,∥ = 1.0 × 10 P 4 / are selected based on the references [35,36]. The parallel current ∥ is calculated by the Ohm's law as the Eq.(6). The curvature drift • (2 × ) and ExB convection *×, • are included in the vorticity equation. The calculated steady state radial electric field (Er) has been compared with experimental measurements in discharge 1080321020 using the pedestal CXRS diagnostic 37 . The simulated Er is shown in Fig. 5 for two cases with and without crossfield drifts. Er is larger in the low field side (LFS) while is smaller in the high field side (HFS) both with and without cross-field drift effects included, which agrees generally with experimental results 38 . Without magnetic drift effects, inside the separatrix, the E ×B flow is balanced by the diamagnetic flow within the pedestal. The Er is mostly determined by ion diamagnetic term ∇R = 5 S $ (F $ ∇ " # in the pedestal with no net flow. While in the SOL region, a positive Er is formed due to the SBCs as shown by the blue curves. Magnetic drift has a significant effect on Er through their contributions to the net flow due to the charge separation 29 . In the simulation of discharge #1080321020, the radial electric field forms a deep (~150kV/m) negative well inside the separatrix due to the steep gradient of density and temperature profiles. The simulation results have similar magnitude and profile shape when compared to experimental measurements shown by the green squares in Fig. 5(b).

The heat flux simulation results
The turbulence (6-field module) simulations under the BOUT++ framework show that in the radial direction, the particle and heat flux are turbulently transported by fluctuations across the separatrix into the SOL, and that parallel transport in the SOL then causes particles and heat to flow from outer midplane towards the divertor targets. From 6-field simulations, we can get the well saturated perpendicular particle and heat flux for both ion and electron at outer midplane. However, due to multiscale nature of boundary plasma turbulent transport, the turbulence simulations are able to reach well saturated electron divertor heat fluxes, but not for ion divertor heat flux 26 . Therefore, we developed a plasma fluid transport model with all drifts included to arrive at a steady state solution for both electrons and ions more efficiently. Fig. 6 shows that both ion and electron parallel heat fluxes reach a steady state at the outer divertor target. The perpendicular heat flux is determined by the cross-field physics governed by drifts and turbulent transport: Where the first terms represent the turbulent flux and the second terms are the drifts flux. (0T" and 1#/" are ExB velocity and diamagnetic velocity, respectively. U ' ∇F is the effective radial velocity from particle diffusion. Fig. 7 shows an example of the perpendicular transport coefficients for the low current case (Ip=0.8 MA). Inside the separatrix, the perpendicular transport coefficients of #" and (" are prescribed in order to match the experimental profiles as described in Eqs. (11)(12) of section 2 and in reference 29 , and the inferred separatrix value is used for the whole SOL region.
Because a core boundary flux Γ #,( | 0G2 = ( "( EH % E0 | 0G2 is unknown, in our simulations, the boundary gradients are calculated from experimental profiles and a nominal thermal diffusivity "#,( =1m 2 /s is given. As shown in table 2, the electron heat transport coefficient near the separatrix is comparable to that calculated from turbulence simulations. In the simulation, the perpendicular heat transport coefficients are the same for ions and electrons ( #" = (" = " ) due to the same ion and electron temperature profiles assumed for C-Mod discharges.  Table 2: Heat transport coefficient at separatrix and heat flux width from transport code and turbulence code for 4 difference current cases.
We impose a sheath boundary condition by assuming that the ion velocity matches the sound speed @ at the divertor targets. The sheath heat flux calculations for both ions and electrons are:  Table 3. The peak of total heat flux qt=qi+qe is close to that measured on C-Mod. The heat flux increases with the current. The electron heat fluxes are dominant at the lower values of Bp (or Ip). However, the ion heat flux become larger with the current increase, approaching to the electron heat flux for the high current case (Ip =1. 4 Table 3: Peak of parallel heat fluxes and widths at outer divertor target for 4 different plasma currents. Here qt is the total parallel heat flux, qi the ion parallel heat flux and qe the electron parallel heat flux. The power crossing the separatrix is the integral of perpendicular heat flux over the separatrix surface and the power at the divertor is the integral of parallel heat flux over the divertor targets. The last two columns in Table 4 show the simulation results of the perpendicular power across the separatrix Qsep and the parallel power at the divertor targets Qtarget. As the current increases, the power across the separatrix increases, which agrees with the experiments as shown in the second column. The simulated power at the divertor targets is also consistent with the simulated power across the separatrix. The value of simulated power at the divertor targets is a little bit smaller than that across separatrix due to energy lost to the chamber wall (<11.5%). The total simulated divertor power is ~2x the experimentally inferred power on the outer divertor. The inner divertor power is not measured for these discharges. Radiative losses in the divertor and SOL are neither measured nor simulated here.  Table 4: The experimental values of Qsol are estimated based on the power balance, which should be considered as upper bounds on the actual Qsol due to the absorption fraction of ICRF being less than unity. The third column is the power deposited on the outer divertor targets from the experiment. Last two columns are the perpendicular power across the separatrix Qsep and the power deposited on the divertor targets Qtarget from the BOUT++ transport simulation.

w/ drifts
In order to compare the BOUT++ transport simulation results of ! with experimental results, we use the same fitting method as that used to populate the experimental heat flux data base 5 , as follows: and ̅ = ( @(& − ) • 0 , where R is the major radius, @(& is the separatrix position at the outer midplane and 0 is the effective flux expansion. The quantity S is the divertor power spreading parameter which depends on local divertor plasma parameters and geometry, ,V is the background heat flux and 2 is the peak heat flux. The fitting function is a convolution from an exponential decay of the outer midplane parallel heat flux profile with a Gaussian function of width S. This formula will be used to calculate the heat flux width ! from experimental data and simulated data. Direct comparisons of simulated heat flux profiles with experimental profiles are shown in Figure 8. Simulations of lower current discharges show a sharp drop in q|| in the private flux region (PFR) as shown by black squares in Fig. 8(a). This does not generally agree with experimental measurements of divertor parallel heat flux as shown by gray squares in Fig. 8(a), which were made with infrared imaging (IR) 39 , which has an estimated spatial resolution in this coordinate system of approximately 1mm. In order to compare the simulated heat flux profile with the measurement, a convolution with an instrument function of 1mm width is applied to the simulated data. Fig. 8(b) shows the overlay of the same IR data with synthetic data points (green squares) by applying the instrument function convolution to the simulation data. This shows a better comparison between the IR measured heat flux profile and synthetic reconstruction based on the BOUT++ simulation. Fig. 8(c) compares a BOUT++ transport simulation and experimental data for the high current case (Ip=1.4 MA). The black squares in Fig.  8(c) are calculated parallel heat flux from BOUT++ simulation, and the light cyan squares and gray squares in Fig. 8(c) are the measured outer divertor parallel heat flux, mapped to the outer midplane. The experimental data come from a suite of divertor diagnostics with high spatial resolution 12 , such as the surface thermocouples 40 for the light cyan squares and Langmuir probes 41 for gray squares. There is the natural scatter that is always present in experimental data. There are clear outliers in the profile and even negative values of heat flux due to both real fluctuations, instrumental reporting error, and uncertainty in the magnetics reconstruction. Here a radial shift of 1mm is applied to bring the peak of the profile close to the strike point R-Rsep=0. The red curve and blue curve in Fig. 8(c) are fits to the profiles using Eq (19). The sharp fall off in q|| for R-Rsep<0 is reproduced well by the simulation, while the simulated heat flux channel for R-Rsep>0 is wider than that in experiment.
The detailed comparison of simulated data with experiment for 4 different currents is listed in table 3. Here we assume that Eq (19) effectively captures the instrumental resolution of the diagnostics in the S parameter, and that the fitted ! can be compared with that fitted to BOUT++ results. Here BOUT++ transport simulation results have an Eich heat flux width ! systematically larger than the Eich scaling would predict. The values are 0.8 to 1.8 times that of the experimental results, which exhibit significant scatter about the projection of the Eich scaling. The electron heat flux width !,( is larger than those of the ions !,# as shown in Fig. 9(a). The same trend is found for both ion and electron heat flux width !,# and !,( , respectively. The transport simulation results for heat flux width is consistent with the turbulence code calculation and follows the Eich scaling trend for low current cases. It should be noted that for the high current case the calculated result does not follow the inverse Ip scaling but is consistent with Goldston HD model as shown by the blue square in Fig. 9(b), possibly due to the high separatrix temperature. Further discussion will be given in the next paragraph. Fig. 8. (a) The overlay of the outer divertor parallel heat flux footprint mapped back to the outer midplane from BOUT++ simulation (black squares) and from IR camera measurement (gray squares); the overlay of the fitting profiles to the simulated data points (red curve) and IR camera data (blue curve) using the Eich-fitting formula with parameters as listed for low current case Ip=0.8 MA; (b) The similar overlay to (a), but BOUT++ simulation data is synthetic reconstruction (green squares) by applying a convolution with an instrument function of 1mm width to the simulation data and the new Eich fitting profile (orange curve) with the parameters as listed; (c) The overlay of the outer divertor parallel heat flux footprint mapped back to the outer midplane from BOUT++ simulation (black squares) and from Langmuir probe(LP) measurement (gray squares) and the surface thermocouples(STC) measurement (light cyan squares); their fitting profiles to the simulated data points (red curve) and experimental data (blue curve) using the Eich-fitting formula with parameters as listed for high current case Ip=1.4 MA; Fig. 9. Heat flux widths for four selected C-Mod discharges with different plasma currents or poloidal magnetic fields Bpol,MP at the outer midplane. (a) the width for ion !,# (red stars), electron !,( (blue stars) and total heat flux !,A (green stars); (b) Simulation data of total heat flux width overlaid on independent fits of experimental heat flux decay length ( ! or =WX ). The black solid curve is the experimental scaling for the ITPA multi-machine database to the Eich fitting formula 9 and black dashed curves for the error bars of the scaling law. Red circles are calculated from 6-field module 26 ; green stars are calculated from the trans-Er module; blue squares are calculated from Goldston's HD model and the gray markers are from the experimental data.
To explore why the simulated heat flux width for the high current case does not continue the trend of Eich scaling, we make a comparison of simulated divertor heat flux width to the Goldston HD model 16 . The model assumes that the magnetic (∇B and curvature B) drifts balance against the near-sonic parallel flow in setting the ion heat flux width, and that cross-field transport at the separatrix is dominated by classical drifts 16,42 .The heuristic drift-based (HD) model gives a width ! ∼ @ , where @ is the ion gyro-radius and q is the safety factor. Goldston's formula for heat flux width follows:  Fig. 9(b) as blue squares. The error bars are calculated from the uncertainty of the separatrix temperature, considering the uncertainty of the separatrix position given by the magnetic equilibrium reconstruction: approximately ~0.01. The short vertical blue lines are the error bars as shown in Fig. 9(b). The simulation results calculated by the transport code under the BOUT++ framework show a similar trend of heat flux width to Goldston's HD model. In the simulation, the averaged parallel flow from outer midplane to divertor target is calculated to be around 0.3-0.4Cs, which is smaller compared with an estimate of 0.5Cs in Goldston's model. Compared to the low current case, the simulated heat flux width becomes larger at the high current. That is because the width ! not only depends on the poloidal magnetic field at the outer midplane Bpol,MP (Bp) but also is very sensitive to the separatrix temperature. In order to estimate the dependence of the separatrix temperature, we do a separatrix temperature scan for high current case (1.4MA). The results are shown in Fig. 10(a). From the transport simulations as shown in Fig. 10(b), it can be found that, the heat flux width increases with the separatrix temperature, which is consistent with the Goldston formula in Eq (16). In Fig. 9(a) the short vertical green line for Bpol,mp=1.1T is the error bar from the scanning of separatrix temperature for high current case as shown in Fig.10(b). When the separatrix temperature is low, the heat flux width is close to the experimental measurement.

The effects of drifts on heat flux width
Classical particle motion from ExB and magnetic drifts is important for understanding tokamak edge/scrape-off-layer (SOL) transport even in the presence of turbulent transport, and this may influence the overall behavior of the heat flux at the divertor targets 42 . The asymmetry of the plasma density and temperature profiles changes between the inner and outer divertor regions and the profiles are poloidally nonuniform with the drift effects [43][44][45] . In this section, the impact of magnetic drift and ExB drift on the divertor heat flux width will be studied by turning on or off the options of magnetic drift and ExB drift. Based on fits to Eq (19), the heat flux width is calculated for the three cases: without any drifts, with the magnetic drift only and with both magnetic and ExB drifts.
The data points of the heat flux width for four C-Mod EDA H-mode discharges with different drifts are shown in Fig. 11. We find that the trend of heat flux width with poloidal field is similar for the cases with or without the drifts. The magnetic drift has a significant influence on the heat flux width. With the magnetic drift, the width is reduced around 2~3 times because the ∇B and curvature B drifts can cause a large parallel flow and hence reduce the SOL residence time. The values calculated are similar to Goldston's HD model, which is also without ExB drift. Reiser and Eich conducted the studies of drift-based scrape-off particle width in X-point geometry to confirm the HD model for the averaged density decay length 42 . In order to know the effect of the ExB drift, the heat flux width is then calculated with both magnetic and ExB drifts as shown by the green stars in Fig. 11. ExB drift can further increase the parallel transport and decreases the heat flux width. From the simulation results, the ExB drift decreases the heat flux width by 10%~25%, which brings the calculated heat flux width closer to the experimental values. Fig. 11. Heat flux width at outer divertor target of four C-Mod EDA H-mode discharges. Blue bullets are the heat flux width without any drifts, red squares are with magnetic drift and green stars are with both magnetic and ExB drifts.

The competition between drifts and turbulent transport for C-Mod
We expect that the heat flux width is determined by the competition between the cross-field and parallel transport 46 . The cross-field transport is governed by drifts and turbulent transport. The effect of drifts with increasing particle turbulent diffusivity has been investigated by reference 47 . In this section, the competition between drifts and turbulent heat transport coefficient is investigated. We perform a scan of the turbulent heat transport coefficients for the 0.8MA C-Mod discharge. The perpendicular heat transport coefficient ( " ) in the scrape-off-layer (SOL) is set constant in radius, and then scanned in a series of simulations. Inside the separatrix, we use the same " interpreted from experimental profiles. The SOL " varies by 160X from 0.05 4 / to 8 4 / , as shown in Fig. 12(a). Fig. 12(b) shows the transport simulation results of heat flux width ! for different " with drifts. The scan of turbulence diffusivity " identifies two distinct regimes. When " < 0.5 4 / , the heat flux width ! is a constant, the drifts dominate cross-field transport and heat flux width is insensitive to the turbulent transport " in this regime, which is consistent with Goldston HD-model as shown by blue star in Fig. 12(b). Goldston HD model yields a lower limit of the width ! . When " > 0.5 4 / , heat flux width ! increases with " , corresponding to turbulence dominating cross-field transport. The " ) = 0.5 4 / is the critical thermal diffusivity value for the transition from drift dominant regime to turbulent dominant regime. The green star shows the interpreted result from the experimental profile using the formula (12) with " = 0.41 4 / and ! = 1.27 mm. The red bullet in Fig. 12(b) is the simulation result from turbulence code with " = 0.23 4 / and ! = 1.04 mm, which is comparable to transport code results. The similar scan of thermal diffusivity has been done for the high current case(1.4MA) and a smaller critical " ) = 0.3 4 / is found due to the decrease of drift effect. Both the thermal diffusivity calculated by turbulence code and that interpreted from the experimental data are smaller than the critical thermal diffusivity 14 . The simulation results suggest in C-Mod EDA H-mode discharges, drifts and turbulent transport compete to determine the characteristic heat flux width, and that drifts set a lower bound on ! . However, uncertainty remains, since the heat flux width calculated by all models systematically overpredict the measured values as shown by the gray dashed dot horizontal line. It should be mentioned that the predictions for future machines, such as CFETR and ITER, show a lower critical thermal diffusivity " ) and a higher width ! .
ITER and CFETR are predicted to be in a turbulence dominant regime 14,15 .

Summary
The transport code (trans-Er) with all the cross-field drifts under the BOUT++ framework is used to calculate the steady state radial electric field (Er) and to study the physics setting heat flux width. Transport coefficients are calculated from the experimental profiles inside the separatrix, then extending to the SOL. Simulation results show that Er is larger in the LFS than that in the HFS. The steady state Er has been compared favorably to experimental measurements from a C-Mod enhanced Dα (EDA) H-mode discharge.
In order to understand the relative role of drifts vs turbulent transport in setting the heat flux width, a set of four C-Mod EDA H-mode discharges with lower single null divertor configuration are simulated. The simulation results yield similar ! to experimental measurements within a factor of 0.8-1.8 and show a trend similar to Goldston's HD model. In simulations, the power across the separatrix increases with current for these four discharges. The power deposited on the divertor target is consistent with the power across the separatrix due to the low SOL radiation loss. Magnetic drift has a significant influence on the heat-flux width, while the ExB drift decreases the heat flux width by 10%~25%, which gives improved agreement with experiment relative to Goldston's model. However, there are indications from modeling a separate EDA H-mode discharge that the simulation with full magnetic drifts may overestimate Er. In simulations of C-Mod discharges, the heat flux width is sensitive to the temperature near the separatrix. A turbulence diffusivity scan ( " ) identifies two distinct regimes: a drift dominant regime when " is small and a turbulence dominant regime when " is large. Magnetic drifts give a lower limit of the width ! , similar to Goldston. This exercise suggests that drifts and turbulent transport compete to determine the heat flux width in C-Mod EDA H-mode discharges.